[ICPC 2019 Nanchang Invitational Online G]tsy’s number

这几周很忙…好久没有静心做点题了。。
北理校赛的期望DP题做不对
南昌网赛的积性函数题不会做
各种专题都有待加强
这几天先再看看数论吧。

Link

Solution

首先

\begin{aligned}
&\frac{\varphi(i)\varphi(j^2)\varphi(k^3)}{\varphi(i)\varphi(j)\varphi(k)}\\
=&\frac {j^2\prod_{p|j^2}(1-\frac 1p)k^3\prod_{p|k^3}(1-\frac 1p)}
{j\prod_{p|j}(1-\frac 1p)k\prod_{p|k}(1-\frac 1p)}\\
=&jk^2
\end{aligned}

式子化为

\begin{aligned}
&\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^njk^2\varphi((i,j,k))\\
=&\sum_{d}\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
(jd)(kd)^2[(i,j,k)=1]\\
=&\sum_{d}d^3\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
jk^2[(i,j,k)=1]\\
=&\sum_{d}d^3\varphi(d)\sum_{i=1}^{[\frac nd]}\sum_{j=1}^{[\frac nd]}\sum_{k=1}^{[\frac nd]}
jk^2\sum_{t|i,t|j,t|k}\mu(t)\\
=&\sum_{d}d^3\varphi(d)\sum_{t}t^3\mu(t)\sum_{i=1}^{[\frac n{dt}]}\sum_{j=1}^{[\frac n{dt}]}\sum_{k=1}^{[\frac n{dt}]}jk^2\\
=&\sum_{T}T^3\sum_{d|T}\varphi(d)\mu(\frac Td)\sum_{i=1}^{[\frac nT]}\sum_{j=1}^{[\frac nT]}\sum_{k=1}^{[\frac nT]}jk^2\\
=&\sum_{T}T^3\sum_{d|T}\varphi(d)\mu(\frac Td)[\frac nT](\sum_{j=1}^{[\frac nT]}j)(\sum_{k=1}^{[\frac nT]}k^2)
\end{aligned}

考虑筛出T^3(\varphi*\mu)(T)的前缀和。
都是积性函数,随便筛的。

\begin{aligned}
&(p^t)^3(\mu * \varphi)(p^t)\\
=&p^{3t}\sum_{d|p^t}\mu(d)\varphi(\frac {p^t}d)\\
=&p^{3t}(\mu(1)\varphi(p^t)+\mu(p)\varphi(p^{t-1}))\\
=&p^{3t}(\varphi(p^t)-\varphi(p^{t-1}))\\
\end{aligned}

Code

线性筛也写的不熟,写一次忘一次。

#include <bits/stdc++.h>

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

unsigned int Sum(unsigned int x) {
    return (long long)x*(x+1)/2;
}

unsigned int Sum2(unsigned int x) {
    return (unsigned long long)x*(x+1)%(6ull<<30)*(2*x+1)%(6ull<<30)/6;
}

unsigned int Pow3(unsigned int x) {
    return x*x*x;
}

unsigned int g[XN];
void Prep() {
    static bool notPrime[XN];
    static struct { int p,pt; }mp[XN];
    static int prime[XN],pcnt;

    g[1]=1;

    for(int i=2;i<=N;++i) {
        if(!notPrime[i]) {
            prime[++pcnt]=i;
            mp[i]={i,i};
            g[i]=Pow3(i)*(i-2);
        }
        for(int j=1,x;j<=pcnt && prime[j]<=mp[i].p && (x=i*prime[j])<=N;++j) {
            notPrime[x]=1;
            if(mp[i].p==prime[j])
                mp[x]={prime[j],mp[i].pt*prime[j]};
            else
                mp[x]={prime[j],prime[j]};
            auto Phi=[&](int pt)->int {
                return pt==1?1:(mp[x].p-1)*(pt/mp[x].p);
            };
            g[x]=g[x/mp[x].pt]*Pow3(mp[x].pt)*(Phi(mp[x].pt)-Phi(mp[x].pt/mp[x].p));
        }
    }
    for(int i=1;i<=N;++i)
        g[i]+=g[i-1];
}

void Work() {
    int n;std::cin>>n;
    unsigned int Ans=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        Ans+=(n/l)*Sum(n/l)*Sum2(n/l)*(g[r]-g[l-1]);
    }
    std::cout<<(Ans%(1<<30))<<'\n';
}

int main() {
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    //freopen("input","r",stdin);
    Prep();
    /*
    for(int i=1;i<=100;++i)
        std::cerr<<Sum(i)-Sum(i-1)<<' '<<Sum2(i)-Sum2(i-1)<<'\n';
    */
    int T;std::cin>>T;
    while(T--)
        Work();
    return 0;
}

[BZOJ 4176]Lucas的数论

Link

Solution

比较套路地画出来

\sum_{d}\mu(d)F([\frac nd])^2

其中

F(n)=\sum_{i=1}^n[\frac ni]

主要是求F函数部分的复杂度如何估计

\begin{aligned}
T(n)\approx&\sum_{d=1}^{\sqrt n}(\sqrt d+\sqrt{\frac nd})\\
\approx&\int_{1}^{\sqrt n}(\sqrt x+\sqrt {\frac nx}){\rm d}x\\
=&(2(nx)^{\frac 12}+\frac {2x^{\frac 32}}3)\bigg|_1^{\sqrt n}\\
=&O(n^{\frac 34})
\end{aligned}

Code

const int P=1000000007,N=1e6,XN=N+11;

int prime[XN],mu[XN],pcnt,smu[XN];
void Prep() {
    static bool notPrime[XN];
    mu[1]=1;
    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];
        }
    }
    for(int i=1;i<=N;++i)
        smu[i]=smu[i-1]+mu[i];
}

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

int F(int n) {
    long long res=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        (res+=(long long)(r-l+1)*(n/l))%=P;
    }
    return res;
}

int main() {
    Prep();
    int n;fin>>n;
    long long Ans=0;
    for(int l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        long long x=F(n/l);
        (x*=x)%=P;
        (Ans+=(long long)(Smu(r)-Smu(l-1))*x)%=P;
    }
    Ans<0 && (Ans+=P);
    fout<<Ans<<'\n';
    return 0;
}

总结一些积性函数的性质

  1. \sum_{d|n}\mu(d)=[n=1]
    x=\prod p_i^{k_i}f(x)=\sum_{d|x}\mu(d)
    f(x)=\prod f(p_i^{k_i})
    单独考虑f(p_i^{k_i})=f(p_i)=\mu(1)+\mu(p_i)=0
    因此f(x)=0
  2. \mu^2(n)=\sum_{d^2|n}\mu(d)
    \mu^2(n)=1\Leftrightarrow n是无平方因子数
    i是无平方因子数时,\sum_{d^2|n}\mu(d)=\mu(1)=1
    否则,\sum_{d^2|n}\mu(d)。现在限定d是无平方因子数。
    n的次数\ge 2的质因子集合为S
    那么,\sum_{d^2|n}\mu(d)=\sum_{P\subset S}(-1)^{|P|}=0
  3. n=\sum_{d|n}\varphi(d)
    不会证明
  4. \sum_{i=1}^n[(n,i)=1]i=\frac {n\varphi(n)+[n=1]}{2}
    如果in互质,则(n,i)=(n,n-i)=1,即与n互质的数字成对出现且和为n。因此,对于n\ge 2的情况,有\frac {\varphi(n)}2对,结论就显然了。
  5. D(i)i的约数个数,则D(ij)=\sum_{x|i}\sum_{y|j}[(x,y)=1]
    建立ij的约数d到有序数对(\gcd(i,d),\frac {d}{\gcd(i,d)})的映射。
    那么,一个与一个因数d构成映射的有序数对(a,b)满足a|i,b|j,\gcd(\frac ia,b)=1。而对于任意满足上述条件的数对,都可以反解出d来。
    因此,d与上述数对构成一一映射。
    \begin{aligned} D(ij)=&\sum_{x|i}\sum_{y|j}[(\frac ix,y)=1]\\ =&\sum_{x|i}\sum_{y|j}[(x,y)=1] \end{aligned}

[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;
}

SPOJ DIVCNTK

Link

Solution

使用min_25筛
由于求h的过程都只会用到[\frac ni]集合中的点,那么求g的过程中把[\frac ni]构成的集合求出来,对于这个集合的所有点分层递推(递推过程中需要用到的g也属于[\frac ni]集合),就可以求出g(i,|P|)
然后递归求h,并不需要记忆化。边界见代码。

Code

typedef unsigned long long qword;

const int N=2e5,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;
        }
    }
}

long long n,k;
int lim,psz;

struct Identifier {
    int id[2][XN*2],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;

qword g[XN*2];

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

void Work() {
    lim=sqrt(n)+0.5;
    id.cnt=0;
    for(int i=0;i<=lim;++i)
        id.id[0][i]=id.id[1][i]=0;
    psz=std::lower_bound(prime+1,prime+1+pcnt,lim)-prime;
    static long long kp[XN*2];int kpc=0;
    for(long long 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 && (long long)prime[j]*prime[j]<=kp[i];++i)
            g[id[kp[i]]]-=g[id[kp[i]/prime[j]]]-(j-1);
    std::cout<<H(n,1)+1<<'\n';
}

int main() {
    Prep();
    int T;fin>>T;
    while(T--) {
        fin>>n>>k;
        Work();
    }
    return 0;
}

胡扯min_25筛

记质数集合为P

Problem

以低于线性的复杂度求一类积性函数的和,函数需满足
f(p),f(p^k)(p\in P)很好求。

Solution

{\rm Min}(x)x的最小质因子。

g

g(n,i)=\sum f(x)[ {x\le n {\ \rm and\ } (x \in P {\ \rm or\ } {\rm Min}(x)>P_i) } ]
n < P_i^2时,g(n,i)=g(n,i-1),因为最小的{\rm Min}(x)= P_i的合数就是P_i^2
否则,考虑从g(n,i-1)中减去{\rm Min}(x)= P_i的合数的贡献,也就是\sum f(x)[{x\le n,{\rm Min}(x)=P_i}],提出f(P_i),式子变成了f(P_i)(g([\frac n{P_i}],i-1)-\sum_{j=1}^{i-1}f(P_j))
\therefore
g(n,i)= \begin{cases} g(n,i-1)&,n < P_i^2\\ g(n,i-1)-f(P_i)(g([\frac n{P_i}],i-1)-\sum_{j=1}^{i-1}f(P_j))&,\text{otherwise} \end{cases}

h

h(n,i)=\sum_{j=1}^nf(j)[{\rm Min}(j)\ge P_i]
对于h(n,i)n以内质数的贡献是g(n,|P|)-\sum_{j=1}^{i-1}f(P_j)
现在考虑合数的贡献。
枚举合数的最小质因子P_j,那么最小质因子为P_j的合数的贡献就是\sum_{t\ge 1,P_j^{t+1}\le n}(h([\frac n{P_j^t}],j+1)f(P_j^t)+f(P_j^{t+1}))
因此
h(n,i)=g(n,|P|)-\sum_{j=1}^{i-1}f(P_j)+\sum_{j\ge i}\sum_{t\ge 1,P_j^{t+1}\le n}(h([\frac n{P_j^t}],j+1)f(P_j^t)+f(P_j^{t+1}))
答案就是h(n,1)了。

Tips

具体实现的时候,g不一定求\sum f(x)[ {x\le n {\ \rm and\ } (x \in P {\ \rm or\ } {\rm Min}(x)>P_i) } ],求出的g只要能使得上式可以被很快求出即可。

[zoj4056][ICPC 2018 Qingdao Online]Press the Button

Link

Solution

问题最后是求一个长度为t的数轴上有多少对相邻的关键点距离\ge v,关键点是所有ac的倍数的点。
假设a < c
如果a < v则不存在
否则
考虑c的倍数与它后面相邻的a的倍数的距离
根据根据鸽巢原理,至多c步,这个距离就会开始循环
然后把t分成三部分:
不在循环+循环部分+一段没循环完的部分
细节巨多写到脑壳疼
~正解是模拟到{\rm lcm}(a,b),直接判断,然后循环~

Code

const int XN=1e6+11;
void Work() {
    long long v,t,a,b,c,d;
    fin>>a>>b>>c>>d>>v>>t;
    long long Ans=(long long)((t/a)+1)*b+(long long)((t/c)+1)*d-1;
    if(a>v && c>v) {
        if(a>c)
            std::swap(a,c);
        static int us[XN],ut;
        static long long crec[XN],prec[XN];
        ++ut;long long cnt=0,pc=0,py=0;
        for(;;pc+=c) {
            long long x=(pc/a)*a,y=x+a; 
            cnt+=((pc-x)>v)+((y-pc)>v)+(x-py)/a;
            py=y;
            if(us[y-pc]==ut)
                break;
            us[y-pc]=ut;
            crec[y-pc]=cnt;
            prec[y-pc]=py;
        }
        if(t>prec[py-pc]) {
            long long circnt=cnt-crec[py-pc],cirlen=py-prec[py-pc];
            Ans-=circnt*((t-prec[py-pc])/cirlen)+crec[py-pc];
            py+=((t-prec[py-pc])/cirlen-1)*cirlen;
            pc=(py/c+1)*c;
        } else {
            py=pc=0;
        }
        if(pc<=t)
            for(;;pc+=c) {
                long long x=(pc/a)*a,y=x+a;
                Ans-=((pc-x)>v)+(x-py)/a;
                if(y<=t) {
                    Ans-=((y-pc)>v);
                }
                py=y;
                if(pc+c>t)
                    break;
            }
        if(py<t)
            Ans-=(t-py)/a;
    }
    fout<<Ans<<'\n';
}

int main() {
    int T;fin>>T;
    while(T--)
        Work();
    return 0;
}