【杜教筛/min25筛】计蒜客 Easy Math

S o u r c e : Source: Source:ACM-ICPC 2018 徐州赛区网络预赛
I d e a : Idea: Idea:
f ( n , m ) = ∑ i = 1 m μ ( i ∗ n ) = μ ( n ) ∗ ∑ i = 1 m μ ( i ) [ g c d ( i , n ) = 1 ] f(n,m)=\sum_{i=1}^{m}\mu(i*n)=\mu(n)*\sum_{i=1}^{m}\mu(i)[gcd(i,n)=1] f(n,m)=i=1mμ(in)=μ(n)i=1mμ(i)[gcd(i,n)=1] = μ ( n ) ∗ ∑ d ∣ n μ ( d ) ∑ i = 1 ⌊ m d ⌋ μ ( i ∗ d ) = μ ( n ) ∗ ∑ d ∣ n μ ( d ) f ( d , ⌊ m d ⌋ ) ( 莫 比 乌 斯 反 演 ) =\mu(n)*\sum_{d|n}\mu(d)\sum_{i=1}^{\lfloor \frac{m}{d} \rfloor} \mu(i*d) = \mu(n)*\sum_{d|n}\mu(d)f(d, \lfloor \frac{m}{d} \rfloor) (莫比乌斯反演) =μ(n)dnμ(d)i=1dmμ(id)=μ(n)dnμ(d)f(d,dm)()
f ( 1 , m ) = ∑ i = 1 m μ ( i ) , f ( n , 1 ) = μ ( n ) ( 杜 教 筛 ) f(1, m) = \sum_{i=1}^m \mu(i), f(n, 1) = \mu(n) (杜教筛) f(1,m)=i=1mμ(i),f(n,1)=μ(n)()
也可以考虑min_25筛,素数部分减去n的质因子,合数部分减去这些质因子的倍数。

C o d e : Code: Code:

//杜教筛
#include<bits/stdc++.h>
using namespace std;
#define Toocold
#define I inline
#define lc o<<1
#define rc o<<1|1
#define fi first
#define se second
#define pb push_back
#define ALL(X) (X).begin(), (X).end()
#define bcnt(X) __builtin_popcountll(X)
#define CLR(A, X) memset(A, X, sizeof(A))
using DB = double;
using LL = long long;
using PII = pair<int, int>;
#ifdef Toocold
#define dbg(args...)\
do { cout << "DEBUG: " << #args << " -> "; err(args); } while(0)
#else
#define dbg(...)
#endif // Toocold
void err() { puts(""); }
template<template<typename...> class T, typename t, typename... Args>
void err(T<t> a, Args... args) { for(auto x:a) cout << x << ' '; err(args...); }
template<typename T, typename... Args>
void err(T a, Args... args) { cout << a << ' '; err(args...); }
/*-----------------------------------------------------------------------------*/

const int N = 1e6+10;

int p[N], pnt;
bool isn[N], vis[N];
LL g[N], mu[N], s[N], n;
I LL getg(LL x) { return x<N?s[x]:g[n/x]; }

void init() {
    isn[1] = 1; mu[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!isn[i]) { p[pnt++] = i; mu[i] = -1; }
        for(int j = 0; j<pnt && i*p[j]<N; j++) {
            isn[i*p[j]] = 1;
            if(i%p[j] == 0) { mu[i*p[j]] = 0; break; }
            mu[i*p[j]] = -mu[i];
        }
    }
    for(int i = 1; i < N; i++) s[i] = mu[i]+s[i-1];
}

void cal(LL x) {
    if(x < N) return;
    int y = n/x; if(vis[y]) return;
    vis[y] = 1; g[y] = 1;
    for(LL l = 2, r; l <= x; l = r+1) {
        LL t = x/l; cal(t), r = x/t;
        g[y] -= (r-l+1)*getg(t);
    }
}

LL f(LL n, LL m, LL u) {
    if(m == 0) return 0;
    if(m == 1) return u;
    if(n == 1) return getg(m);
    LL ret = 0;
    for(LL d = 1; d*d <= n; d++) if(n%d == 0) {
        LL nu = mu[d];
        ret += nu*f(d, m/d, nu);
        nu = u/mu[d];
        ret += nu*f(n/d, m/(n/d), nu);
    }
    return ret*u;
}

I void work() {
    LL m, n, tn;
    scanf("%lld%lld", &m, &n);
    ::n = m; tn = n;
    cal(m);

    int u = 1;
    if(n < N) u = mu[n];
    else {
        for(int i = 0; 1LL*p[i]*p[i] <= n; i++) if(n%p[i] == 0) {
            n /= p[i]; u *= -1;
            if(n%p[i] == 0) {
                u = 0;
                break;
            }
        }
        if(n != 1) u *= -1;
    }
    if(u == 0) puts("0");
    else printf("%lld\n", f(tn, m, u));
}

int main() {
    init();
    work();
    return 0;
}

//min_25
#include<bits/stdc++.h>
using namespace std;
//#define Toocold
#define I inline
#define lc o<<1
#define rc o<<1|1
#define fi first
#define se second
#define pb push_back
#define ALL(X) (X).begin(), (X).end()
#define bcnt(X) __builtin_popcountll(X)
#define CLR(A, X) memset(A, X, sizeof(A))
using DB = double;
using LL = long long;
using PII = pair<int, int>;
#ifdef Toocold
#define dbg(args...)\
do { cout << "DEBUG: " << #args << " -> "; err(args); } while(0)
#else
#define dbg(...)
#endif // Toocold
void err() { puts(""); }
template<template<typename...> class T, typename t, typename... Args>
void err(T<t> a, Args... args) { for(auto x:a) cout << x << ' '; err(args...); }
template<typename T, typename... Args>
void err(T a, Args... args) { cout << a << ' '; err(args...); }
/*-----------------------------------------------------------------------------*/

const int N = 1e6+10;

LL tmp;
int p[N], pnt;
bool isn[N];
vector<LL> fac;

void init() {
    isn[1] = 1;
    for(int i = 2; i < N; i++) {
        if(!isn[i]) p[pnt++] = i;
        for(int j = 0; j<pnt && i*p[j]<N; j++) {
            isn[i*p[j]] = 1;
            if(i%p[j] == 0) break;
        }
    }
}

namespace min_25 {
    bool isn[N];
    int B, pnt, cnt, idl[N], idr[N];
    LL p[N], val[N], F[N], n;
    LL H[N], pH[N], G[N], pG[N];

    I bool find(const LL &x) {
        return binary_search(ALL(fac), x);
    }

    I LL f(const LL &x, const LL &k) {
        if(k == 1) return find(x)?0:-1;
        return 0;
    } //f(x^k) = ?

    I void init(LL _n) {
        n = _n; B = sqrt(n+0.5); pnt = cnt = 0;
        for(int i = 2; i <= B; i++) {
            if(!isn[i]) {
                p[++pnt] = i;
                F[pnt] = F[pnt-1]+f(i, 1);
                pG[pnt] = pG[pnt-1]+1;
            }
            for(int j = 1; j<=pnt && p[j]*i<=B; j++) {
                isn[i*p[j]] = 1;
                if(i%p[j] == 0) break;
            }
        }
        for(LL l = 1, r; l <= n; l = r+1) {
            LL x = n/l; r = n/x;
            if(x <= B) idl[x] = ++cnt; else idr[r] = ++cnt;
            val[cnt] = x;
            G[cnt] = x-1;
        }
        for(int i = 1; i <= pnt; i++)
        for(int j = 1; j <= cnt; j++) {
            if(p[i]*p[i] > val[j]) break;
            LL now = val[j]/p[i];
            now = now>B?idr[n/now]:idl[now];
            G[j] = G[j]-(G[now]-pG[i-1]);
        }
        for(int i = 1; i <= cnt; i++) {
            G[i] = -(G[i]-(upper_bound(ALL(fac), val[i])-fac.begin()));
            dbg(G[i]);
        }
    }

    LL S(LL x, int k) {
        if(x<=1 || p[k]>x) return 0;
        LL ans = (G[x<=B?idl[x]:idr[n/x]]-F[k-1]);
        for(int i = k; i<=pnt && p[i]*p[i]<=x; i++) if(!find(p[i]))
        for(LL t = p[i], e = 1; t*p[i] <= x; t*=p[i], e++)
            ans = (ans+S(x/t, i+1)*f(p[i], e)+f(p[i], e+1));
        return ans;
    }

    I LL sum(LL x) {
        init(x);
        return S(x, 1)+1;
    }
}

I void work() {
    LL m, n; scanf("%lld%lld", &m, &n);
    int u = 1;
    for(int i = 0; 1LL*p[i]*p[i] <= n; i++) if(n%p[i] == 0) {
        n /= p[i]; u *= -1;
        fac.pb(p[i]);
        if(n%p[i] == 0) { u = 0; break; }
    }
    if(n != 1) u *= -1, fac.pb(n);
    printf("%lld\n", u*min_25::sum(m));
}

int main() {
    init();
    work();
    return 0;
}


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值