【LOJ#572】Misaka Network 与求和（莫比乌斯反演，杜教筛，min_25筛）

【LOJ#572】Misaka Network 与求和（莫比乌斯反演，杜教筛，min_25筛）

题面

LOJ
\[ans=\sum_{i=1}^n\sum_{j=1}^n f(gcd(i,j))^k\]
其中\(f(x)\)表示\(x\)的次大质因子。

题解

这个数据范围不是杜教筛就是\(min\_25\)筛了吧。。。
看到次大质因子显然要\(min\_25\)筛了吧。。。
莫比乌斯反演的部分比较简单，懒得写过程了。
\[ans=\sum_{T=1}^n [\frac{n}{T}]^2\sum_{d|T}f(d)^k \mu(\frac{T}{d})\]
后面带个指数好麻烦啊，就假装\(f(x)=f(x)^k\)吧。。。
显然要求的就是\(f\)和\(\mu\)狄利克雷卷积的前缀和。。。
令\(\displaystyle S(n)=\sum_{i=1}^n (f*\mu)(i)\)，一脸杜教筛的感觉，类似杜教筛来写式子。
\[g(1)S(n)=\sum_{i=1}^n (f*\mu*g)(i)-\sum_{i=2}^n g(i)S([\frac{n}{i}])\]
看到\(\mu\)了，直接令\(g(x)=1\)，\((f*u*g)(i)=(f*(u*1))(i)=(f*e)(i)=f(i)\)。
写出来就是：
\[S(n)=\sum_{i=1}^nf(i)-\sum_{i=2}^n S([\frac{n}{i}])\]
然后考虑怎么求\(\displaystyle \sum_{i=1}^n f(i)\)，一脸\(min\_25\)筛。
行，本来以为不是\(min\_25\)筛就是杜教筛，没想到两个一起来。
好了，实现啥的就可以看看代码了。
复杂度因为杜教筛不能提前筛好一部分前缀和，所以似乎是\(O(n^{3/4})\)？？？
不太会算复杂度，那就当做\(O(\mbox{跑得过})\)了。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<map>
using namespace std;
#define ll long long
#define uint unsigned int
#define MAX 100000
int K,blk;uint n;
int w[MAX],id1[MAX],id2[MAX],m;
int pri[MAX],tot;
bool zs[MAX];
uint g[MAX],prik[MAX];
uint fpow(uint a,int b)
{
    uint s=1;
    while(b){if(b&1)s*=a;a*=a;b>>=1;}
    return s;
}
int getid(int x){return (x<=blk)?id1[x]:id2[n/x];}
void pre(int n)
{
    for(int i=2;i<=n;++i)
    {
        if(!zs[i])pri[++tot]=i,prik[tot]=fpow(i,K);;
        for(int j=1;j<=tot&&i*pri[j]<=n;++j)
        {
            zs[i*pri[j]]=true;
            if(i%pri[j]==0)break;
        }
    }
}
uint calc(int x,int y)
{
    if(x<=1||pri[y]>x)return 0;
    uint ret=(g[getid(x)]-y+1)*prik[y-1];
    for(int i=y;i<=tot&&1ll*pri[i]*pri[i]<=x;++i)
    {
        ll t1=pri[i],t2=1ll*pri[i]*pri[i];
        for(int e=1;t2<=x;++e,t1=t2,t2*=pri[i])
            ret+=calc(x/t1,i+1)+prik[i];
    }
    return ret;
}
uint M[MAX];bool vis[MAX];
uint S(int n)
{
    if(vis[getid(n)])return M[getid(n)];
    uint ret=calc(n,1)+g[getid(n)];
    for(int i=2,j;i<=n;i=j+1)
        j=n/(n/i),ret-=(j-i+1)*S(n/i);
    vis[getid(n)]=true;
    return M[getid(n)]=ret;
}
int main()
{
    scanf("%u%d",&n,&K);pre(blk=sqrt(n));
    for(uint i=1,j;i<=n;i=j+1)
    {
        j=n/(n/i);w[++m]=n/i;g[m]=w[m]-1;
        if(w[m]<=blk)id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }
    for(int j=1;j<=tot&&1ll*pri[j]*pri[j]<=n;++j)
        for(int i=1;i<=m&&1ll*pri[j]*pri[j]<=w[i];++i)
            g[i]-=g[getid(w[i]/pri[j])]-(j-1);
    uint ans=0,lt=0,nw;
    for(uint i=1,j;i<=n;i=j+1)
    {
        j=n/(n/i);nw=S(j);
        ans+=(uint)1*(n/i)*(n/i)*(nw-lt);
        lt=nw;
    }
    printf("%u\n",ans);
    return 0;
}

转载于:https://www.cnblogs.com/cjyyb/p/10170630.html