[ICPC 2018 Shenyang Online]Convex Hull

Link

Solution

(n+1)\sum_{i=1}^n\mu^2(i)i^2-\sum_{i=1}^n\mu^2(i)i^3

Min25筛过不去……
根据\mu^2(n)=\sum_{d^2|n}\mu(d)

\begin{aligned}&\sum_{i=1}^ni^2\mu^2(i)\\=&\sum_{i=1}^ni^2\sum_{d^2|i}\mu(d)
\\=&\sum_{d}\mu(d)d^4\sum_{i=1}^{[\frac n{d^2}]}i^2\end{aligned}
\begin{aligned}&\sum_{i=1}^ni^3\mu^2(i)\\=&\sum_{i=1}^ni^3\sum_{d^2|i}\mu(d)
\\=&\sum_{d}\mu(d)d^6\sum_{i=1}^{[\frac n{d^2}]}i^3\end{aligned}

Code

const int N=1e5,XN=N+11;

long long b20=(1ll<<20)-1,P;

long long Add(long long x,long long const &y) {
    return (x+=y)>=P?x-P:x;
}

void AddBy(long long &x,long long const &y) {
    (x+=y)>=P && (x-=P);
}

long long Minus(long long x,long long const &y) {
    return (x-=y)<0?x+P:x;
}

void MinusBy(long long &x,long long const &y) {
    (x-=y)<0 && (x+=P);
}

long long Mul(long long const &x,long long const &y,long long P=::P) {
    return ((((((x>>20)*(y>>20)<<20)+(x>>20)*(y&b20)+(y>>20)*(x&b20))%P)<<20)+(x&b20)*(y&b20))%P;
}

void MulBy(long long &x,long long const &y) {
    x=Mul(x,y);
}

int prime[XN],pcnt,mu[XN];
void Prep() {
    mu[1]=1;
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            } else
                mu[i*prime[j]]=-mu[i];
        }
    }
}

long long Sum2(long long x) {
    return Mul(Mul(x,x+1,P*6),x*2+1,P*6)/6;
}

long long Sum3(long long x) {
    x=Mul(x,x+1,P*2)/2;
    return Mul(x,x);
}

int main() {
    Prep();
    long long n;
    while(std::cin>>n>>P) {
        long long Ans,a=0,b=0;
        for(long long d=1;d*d<=n;++d)
            if(mu[d]) {
                long long sqd=d*d%P,tbd=sqd*d%P;
                AddBy(a,Mul(mu[d]==1?mu[d]:P-1,Mul(Mul(sqd,sqd),Sum2(n/(d*d)))));
                AddBy(b,Mul(mu[d]==1?mu[d]:P-1,Mul(Mul(tbd,tbd),Sum3(n/(d*d)))));
            }
        Ans=Minus(Mul(n+1,a),b);
        fout<<Ans<<'\n';
    }
    return 0;
}

Code(TLE)

const int N=1e5,XN=N+11;

long long b20=(1ll<<20)-1,P;

long long Add(long long x,long long const &y) {
    return (x+=y)>=P?x-P:x;
}

void AddBy(long long &x,long long const &y) {
    (x+=y)>=P && (x-=P);
}

long long Minus(long long x,long long const &y) {
    return (x-=y)<0?x+P:x;
}

void MinusBy(long long &x,long long const &y) {
    (x-=y)<0 && (x+=P);
}

long long Mul(long long const &x,long long const &y,long long P=::P) {
    return ((((((x>>20)*(y>>20)<<20)+(x>>20)*(y&b20)+(y>>20)*(x&b20))%P)<<20)+(x&b20)*(y&b20))%P;
}

void MulBy(long long &x,long long const &y) {
    x=Mul(x,y);
}

int prime[XN*2],pcnt;
void Prep() {
    static bool notPrime[XN*2];
    for(int i=2;i<=N*2;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N*2;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

namespace Min25 {
    std::function<long long(int,int)> F;

    long long n;
    int lim,psz;

    struct Identifier {
        int id[2][XN],cnt;
        int &operator [](long long x) {
            int &res=x<=lim?id[0][x]:id[1][n/x];
            if(res==0)
                res=++cnt;
            return res;
        }
    }id;

    long long g[XN*2],fps[XN];

    long long H(long long n,int m) {
        if(n<=1 || m>psz)
            return 0;
        long long res=Minus(g[id[n]],fps[m-1]);
        for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
            long long pt=prime[i],pt1=pt*prime[i];
            for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
                AddBy(res,Add(Mul(F(prime[i],t),H(n/pt,i+1)),F(prime[i],t+1)));
        }
        return res;
    }

    long long Solve(long long n,std::function<long long(int,int)> F,std::function<long long(long long)> gInit) {
        static long long kp[XN*2];
        int kpc=0;
        lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
        for(int i=id.cnt=0;i<=lim;++i)
            id.id[0][i]=id.id[1][i]=0;
        Min25::F=F;
        Min25::n=n;
        for(long long l=1,r;l<=n;l=r+1) {
            r=n/(n/l);
            g[id[kp[++kpc]=n/l]]=gInit(n/l);
        }
        for(int i=1;i<=psz;++i)
            fps[i]=Add(fps[i-1],F(prime[i],1));
        for(int j=1;j<=psz;++j)
            for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i)
                MinusBy(g[id[kp[i]]],Mul(F(prime[j],1),Minus(g[id[kp[i]/prime[j]]],fps[j-1])));//XXX
        return H(n,1);
    }
}

long long Sum2(long long x) {
    return Mul(Mul(x,x+1,P*6),x*2+1,P*6)/6;
}

long long Sum3(long long x) {
    x=Mul(x,x+1,P*2)/2;
    return Mul(x,x);
}

int main() {
    Prep();
    long long n;
    while(std::cin>>n>>P) {
        //(n+1)*mu^2(i)*i^2-mu^2(i)*i^3
        long long a=Add(1,Min25::Solve(n,[](int p,int t)->long long { return t==1?Mul(p,p):0; },
                                         [](long long x)->long long { return Minus(Sum2(x),1);})),
                  b=Add(1,Min25::Solve(n,[](int p,int t)->long long { return t==1?Mul(p,Mul(p,p)):0; },
                                         [](long long x)->long long { return Minus(Sum3(x),1);})),
                  Ans=Minus(Mul(n+1,a),b);
        fout<<Ans<<'\n';
    }
    return 0;
}

[loj 6053]简单的函数

Link

Solution

p xor 1并不是积性函数
但是除了p=2的时候,p^1=p-1
p和1都是积性函数,就可以筛出来然后加上2把p=2的情况修正好。

Code

const int P=1e9+7,inv2=500000004;

int Add(int x,int const &y) {
    return (x+=y)>=P?x-P:x;
}

int Minus(int x,int const &y) {
    return (x-=y)<0?x+P:x;
}

int Mul(long long x,int const &y) {
    return x*y%P;
}


const int N=2e5,XN=N+11;

int prime[XN],pcnt,psum[XN],fsum[XN];
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
    for(int i=1;i<=pcnt;++i)
        psum[i]=Add(psum[i-1],prime[i]);
}

int lim,psz;long long n;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](long long x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;


int F(int p,int t) {
    return p^t;
}

int g[XN*2];

int H(long long n,int m) {
    if(n<=1 || m>psz)
        return 0;
    int res=Minus(g[id[n]],fsum[m-1]);
    for(int i=m;i<=psz && (long long)prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res=Add(res,Add(Mul(F(prime[i],t),H(n/pt,i+1)),F(prime[i],t+1)));
    }
    return res;
}

int main() {
    Prep();
    static long long kp[XN*2];
    int kpc=0;
    fin>>n;
    lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
    static int g1[XN*2],g2[XN*2];
    for(long long l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        kp[++kpc]=n/l;
        g1[id[n/l]]=Minus((n/l%P)*(n/l%P+1)%P*inv2%P,1);
        g2[id[n/l]]=Minus(n/l%P,1);
    }
    for(int j=1;j<=psz;++j)
        for(int i=1;i<=kpc && (long long)prime[j]*prime[j]<=kp[i];++i) {
            g1[id[kp[i]]]=Minus(g1[id[kp[i]]],Mul(prime[j],Minus(g1[id[kp[i]/prime[j]]],psum[j-1])));
            g2[id[kp[i]]]=Minus(g2[id[kp[i]]],Minus(g2[id[kp[i]/prime[j]]],j-1));
        }
    for(int i=1;i<=psz;++i)
        fsum[i]=Add(fsum[i-1],prime[i]^1);
    for(int i=1;i<=kpc;++i)
        g[i]=Add(Minus(g1[i],g2[i]),kp[i]>=2?2:0);
    fout<<Add(H(n,1),1)<<'\n';
    return 0;
}

[ICPC 2018 Xuzhou Online]Easy Math

Link

Min25

Solution

\sum_{i=1}^n\mu(ix)=\mu(x)\sum_{i=1}^n\mu(i)[\gcd(i,x)=1]

考虑计算

\begin{aligned}
&\sum_{i=1}^nf(i)\\
&f(i)=\mu(i)[\gcd(i,n)=1]
\end{aligned}

首先,

f(i)=\begin{cases}
1,&i=1\\
-1,&i=p{\ \rm and\ }p\not|x\\
0,&i=p^k{\ \rm or\ }p|x
\end{cases}

然后,分\gcd(i,x)的不同情况讨论,容易知道f(i)是积性函数。
因此可以使用min25筛法了。

Tips

还是不必拘泥于g的形式……

Code

const int N=5e4,XN=N+11;

int prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    for(int i=2;i<=N;++i) {
        if(!notPrime[i])
            prime[++pcnt]=i;
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0)
                break;
        }
    }
}

int n,lim,psz;long long x;
struct Identifier {
    int id[2][XN],cnt;
    int &operator [](int x) {
        int &res=x<=lim?id[0][x]:id[1][n/x];
        if(res==0)
            res=++cnt;
        return res;
    }
}id;

long long dv[100];int dt;
int Mu(long long x) {
    int res=1;long long cx=x;
    for(long long d=2;d*d<=cx && x!=1;++d)
        if(x%d==0) {
            x/=d;
            dv[++dt]=d;
            if(x%d==0)
                return 0;
            else
                res*=-1;
        }
    if(x!=1) {
        res*=-1;
        dv[++dt]=x;
    }
    return res;
}

int F(int p,int t) {
    return t==1 && x%p?-1:0;
}

int g[XN*2],fsum[XN];

long long H(int n,int m) {
    if(n<=1 || m>psz)
        return 0;
    long long res=g[id[n]]-fsum[m-1];
    for(int i=m;i<=psz && prime[i]*prime[i]<=n;++i) {
        long long pt=prime[i],pt1=pt*prime[i];
        for(int t=1;pt1<=n;++t,pt=pt1,pt1*=prime[i])
            res+=F(prime[i],t)*H(n/pt,i+1)+F(prime[i],t+1);
    }
    return res;
}

int main() {
    Prep();
    static int kp[XN*2],kpc;
    fin>>n>>x;
    int mun=Mu(x);
    if(mun==0)
        fout<<0<<'\n';
    else {
        lim=sqrt(n)+0.5,psz=std::upper_bound(prime+1,prime+1+pcnt,lim)-prime;
        for(int l=1,r;l<=n;l=r+1) {
            r=n/(n/l);
            g[id[kp[++kpc]=n/l]]=n/l-1;
        }
        for(int j=1;j<=psz;++j)
            for(int i=1;i<=kpc && prime[j]*prime[j]<=kp[i];++i)
                g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);//(1-(j-Count(j)));
        for(int i=1;i<=psz;++i)
            fsum[i]=fsum[i-1]+F(prime[i],1);
        for(int i=1;i<=kpc;++i)
            g[i]=std::upper_bound(dv+1,dv+1+dt,kp[i])-dv-1-g[i];
        fout<<mun*(H(n,1)+1)<<'\n';
    }
    return 0;
}

递归求解+杜教筛

网赛的时候写出了这个式子,但是感觉不大能过就一直在自闭……

Solution

\begin{aligned}
&\sum_{i=1}^n\mu(im)\\=&\mu(m)\sum_{i=1}^n\mu(i)[(i,m)=1]\\=&\mu(m)\sum_{i=1}^n\mu(i)\sum_{d|i,d|m}\mu(d)\\=&\mu(m)\sum_{d|m}\mu(d)\sum_{i=1}^{[\frac nd]}\mu(id)\end{aligned}

S(n,m)=\sum_{d|m}\mu(d)S([\frac nd],d)

Tips

  1. 搞清楚哪个是要算的函数……一开始多乘了一个\mu愣是没看出来。
  2. 这种递归的式子,与其对着自闭,不如大胆写,说不定就过了呢。

Code

const int N=2e6,XN=N+11;

int mert[XN],prime[XN],pcnt;
void Prep() {
    static bool notPrime[XN];
    mert[1]=1;
    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mert[i]=-1;
        }
        for(int j=1;j<=pcnt && i*prime[j]<=N;++j) {
            notPrime[i*prime[j]]=1;
            if(i%prime[j]==0) {
                mert[i*prime[j]]=0;
                break;
            } else
                mert[i*prime[j]]=-mert[i];
        }
    }
    for(int i=2;i<=N;++i)
        mert[i]+=mert[i-1];
}

long long Mert(int n) {
    if(n<=N)
        return mert[n];
    else {
        static std::unordered_map<int,long long> map;
        if(map.count(n))
            return map[n];
        long long &res=map[n]=1;
        for(int l=2,r;l<=n;l=r+1) {
            r=n/(n/l);
            res-=Mert(n/l)*(r-l+1);
        }
        return res;
    }
}

int Mu(long long x) {
    if(x<=N)
        return mert[x]-mert[x-1];
    else {
        static std::unordered_map<long long,int> map;
        if(map.count(x))
            return map[x];
        int &mu=map[x]=1;
        long long m=x;
        for(int i=1;(long long)prime[i]*prime[i]<=m && x!=1;++i)
            if(x%prime[i]==0 && (mu=-mu,(x/=prime[i])%prime[i]==0))
                return mu=0;
        if(x!=1)
            mu=-mu;
        return mu;
    }
}

long long S(int n,long long m) {
    if(n==0)
        return 0;
    else if(m==1)
        return Mert(n);
    else {
        static std::map<std::pair<int,long long>,long long> map;
        if(map.count(std::pair<long long,int>(n,m)))
            return map[std::pair<long long,int>(n,m)];
        long long &res=map[std::pair<long long,int>(n,m)];
        if(Mu(m)==0)
            return res=0;
        for(long long d=1;d*d<=m;++d)
            if(m%d==0) {
                res+=Mu(d)*S(n/d,d);
                if(m!=d*d)
                    res+=Mu(m/d)*S(n/(m/d),m/d);
            }
        return res*=Mu(m);
    }
}

int main() {
    Prep();
    long long m;int n;
    fin>>n>>m;
    fout<<S(n,m)<<'\n';
    return 0;
}