标签：P7504 frac R1 int sum return mu prim HMOI

好像被我推复杂了...

简述题意

求：

\[\sum^{n}_{x=1}\sum^{n}_{y=1}\sum^{x}_{i=1}i[(x,k_{1})=1][(i,x)=1]\sum^{y}_{j=1}[(y,k_{2})=1][(y,j)=1]j \]

\(n,k_1,k_2\leq 2\times 10^9\)。

Solution

首先是推式子：

\[\begin{aligned} ans&=\sum^{n}_{x=1}\sum^{x}_{i=1}i[(x,k_{1})=1][(i,x)=1]\sum^{n}_{y=1}\sum^{y}_{j=1}[(y,k_{2})=1][(y,j)=1]j\\ \end{aligned} \]

注意到左右两边形式一样，设 \(S(n,k)=\sum^{n}_{x=1}\sum^{x}_{i=1}i[(x,k)=1][(i,x)=1]\)。

那么有：

\[\begin{aligned} ans&=S(n,k_{1})\times S(n,k_{2})\\ \\ S(n,k)&=\sum^{n}_{x=1}\sum^{x}_{i=1}i[(x,k)=1][(i,x)=1]\\ &=\sum^{n}_{x=1}[(x,k)=1]\sum^{x}_{i=1}i\sum_{g|\gcd(i,x)}\mu(g)\\ &=\sum^{n}_{x=1}[(x,k)=1]\sum^{}_{g|x}\mu(g)g C_{2}(\frac{x}{g})\\ &=\sum^{n}_{x=1}\sum_{t|\gcd(x,k_1)}\mu(t)\sum^{}_{g|x}\mu(g)g C_{2}(\frac{x}{g})\\ &=\frac{1}{2}\sum^{}_{t|k_{1}}\mu(t)\sum^{\lfloor\frac{n}{t}\rfloor}_{x=1}x\sum^{}_{g|x}\mu(g)\frac{x}{g}+\mu(g)\\ &=\frac{1}{2}\sum^{}_{t|k_{1}}\mu(t)\sum^{\lfloor\frac{n}{t}\rfloor}_{x=1}x\varphi(x)+[x=1]\\ &=1+\frac{1}{2}\sum^{}_{t|k_{1}}\mu(t)\sum^{\lfloor\frac{n}{t}\rfloor}_{x=1}tx\varphi(tx))\\ \end{aligned} \]

如果是求 \(i\varphi(i)\) 的前缀和，可以直接上 Min_25 ，可是这里的形式比较诡异，不能直接求。

不过注意到 \(\mu(t)\not = 0\)，当且仅当 \(t\) 中的每一个质因子只出现一次，似乎是可以做一个DP。

我们令 \(p\) 表示 \(t\) 中的某一个质因子，同时定义一个新的函数 \(G(n,t)\)。

\[\begin{aligned} G(n,t)&=\sum^{n}_{x=1}x\varphi(tx)\\ &=\sum^{n}_{x=1}[\gcd(x,p)=1]x\varphi(x\frac{t}{p})(p-1)+\sum^{\lfloor\frac{n}{p}\rfloor}_{x=1}xp^2\varphi(xp\frac{t}{p})\\ &=\sum^{n}_{x=1}x\varphi(x\frac{t}{p})(p-1)+\sum^{\lfloor\frac{n}{p}\rfloor}_{x=1}xp\varphi(xt)\\ &=(p-1)G(n,\frac{t}{p})+pG(\lfloor\frac{n}{p}\rfloor,t)\\ \\ G(n,1)&=\sum^{n}_{i=1}i\varphi(i)=F(n)\\ \\ S(n,k)&=1+\frac{1}{2}\sum_{t|k \And \mu(t)\not = 0}t\mu(t)G(\lfloor \frac{n}{t} \rfloor,t) \end{aligned} \]

加上记忆化之后，复杂度便正确了。

首先是 \(G\) 部分，总共状态个数是 \(O(k^{\frac{1}{3}}n^{\frac{1}{2}})\) 的，每次转移是 \(O(1)\)。

然后是 \(F\) 部分，只有 \(O(\sqrt{n})\) 个取值，只需要在 Min_25 内加一个记忆化也是可以做到 $O(n^{\frac{3}{4}}) $ 或者 \(O(\frac{n^{\frac{3}{4}}}{\ln n})\) 。

视 \(n,k\) 同阶，那么总复杂度为 \(O(n^{\frac{5}{6}})\) ，实测跑的还是比较快的，完全跑不满，瓶颈反而在 Min_25 内记忆化所使用的 map 中。

Code

#include<cstdio>
#include<iostream>
#include<vector>
#include<algorithm>
#include<cstring>
#include<cmath>
#include<map>
#include<unordered_map>
using namespace std;
#define ll long long 
#define ri int
#define pii pair<int,int>
const ll mod=998244353,inv=mod+1>>1,inv3=332748118;
ll add(ll x,ll y){return (x+=y)<mod?x:x-mod;}
ll dec(ll x,ll y){return (x-=y)<0?x+mod:x;}
ll ksm(ll d,ll t){ll res=1;for(;t;t>>=1,d=d*d%mod) if(t&1) res=res*d%mod;return res;}
const int MAXN=1e5+7;
int totcnt=0;
namespace Getphi{
    const int MAXM=1e7+7;
    int val[MAXM],as[MAXM],nxt[MAXM],mapcnt;
    struct hmap{
        static const int P=10007;
        int hed[P];
        bool count(int x){
            int c=hed[x%P];
            while(c){
                if(val[c]==x) return 1;
                c=nxt[c];
            }return 0;
        }
        int& operator [](int x){
            int c=hed[x%P];
            while(c){
                if(val[c]==x) return as[c];
                c=nxt[c];
            }
            ++mapcnt;val[mapcnt]=x;nxt[mapcnt]=hed[x%P];hed[x%P]=mapcnt;
            return as[mapcnt];
        }
        int size(){return mapcnt;}
    }S[4650];
    int n,p1[MAXN],p2[MAXN],cnt,N,flag[MAXN],prim[MAXN],w[MAXN];
    int pos(int x){
        if(1ll*x*x<=n) return p1[x];
        else return p2[n/x];
    }
    ll g1[MAXN],g2[MAXN];
    ll C2(ll x){return x*(x+1)/2%mod;}
    ll C3(ll x){return C2(x)*(x+x+1)%mod*inv3%mod;}
    void init(int _n){
        n=_n;
        N=sqrt(n);
        for(ri i=2;i<=N;++i){
            if(!flag[i]) prim[++prim[0]]=i;
            for(ri j=1;j<=prim[0]&&i*prim[j]<=N;++j){
                flag[i*prim[j]]=1;
                if(i%prim[j]==0) break;
            }
        }
        for(ri l=1,r;l<=n;l=r+1){
            int cur=n/l;r=n/cur;
            w[++cnt]=cur;
            if(1ll*cur*cur<=n) p1[cur]=cnt;
            else p2[n/cur]=cnt;
            g1[cnt]=C2(cur)-1;
            g2[cnt]=C3(cur)-1;
        }
        for(ri t=1;t<=prim[0];++t){
            for(ri i=1;i<=cnt&&w[i]>=prim[t]*prim[t];++i){
                int nxtl=pos(w[i]/prim[t]),nxtr=(t==1?0:pos(prim[t-1]));
                g1[i]=dec(g1[i],dec(g1[nxtl],g1[nxtr])*prim[t]%mod);
                g2[i]=dec(g2[i],dec(g2[nxtl],g2[nxtr])*prim[t]%mod*prim[t]%mod);
            }
        }
    }
    ll getS(int n,int cur){
        int P=pos(n);
        if(cur&&prim[cur]>=n) return 0;
        if(S[cur].count(P)) return S[cur][P];
        int r=pos(n),l=cur==0?0:pos(prim[cur]);
        ll ans=dec(dec(g2[r],g2[l]),dec(g1[r],g1[l]));
        for(ri i=cur+1;i<=prim[0]&&1ll*prim[i]*prim[i]<=n;++i){
            for(ll e=1,v=prim[i],W=1;v<=n;++e,v*=prim[i],W*=prim[i]){
                ll res=getS(n/v,i);
                if(e!=1) ++res;
                res=res*v%mod*W%mod*(prim[i]-1)%mod;
                ans=add(ans,res);
            }
        }
        return S[cur][P]=ans;
    }
    ll getphi(ll x){
        ll now=getS(x,0)+1;
        return now;
    }
}
ll ans;
ll sta[64],top;
int f[64][MAXN];
ll g(ll n,int cur){
    if(n==0) return 0;
    int P=Getphi::pos(n);
    if(~f[cur][P]) return f[cur][P];
    if(cur==0) return f[cur][P]=Getphi::getphi(n);
    return f[cur][P]=add((sta[cur]-1)*g(n,cur-1)%mod,sta[cur]*g(n/sta[cur],cur)%mod);
}
ll S(ll n,ll k){
    ll res=1;
    for(ri t=1;t*t<=k;++t){
        if(k%t) continue;
        if(t<=n){
            ll x=t;
            int now=1;
            top=0;
            for(ll i=2;i*i<=t&&i<=x;++i){
                int cnt=0;
                while(x%i==0) ++cnt,x/=i;
                if(cnt>1) goto Nex;
                if(cnt) sta[++top]=i,now=-now;
            }
            if(x>1) now=-now,sta[++top]=x;
            for(ri i=0;i<=top;++i) for(ri j=0;j<=Getphi::cnt;++j) f[i][j]=-1;
            if(now== 1) res=add(res,g(n/t,top)*t%mod);
            if(now==-1) res=dec(res,g(n/t,top)*t%mod);
            Nex:;
        }
        if(t*t!=k&&k/t<=n){
            ll x=(k/t);
            int now=1;
            top=0;
            for(ll i=2;i*i<=(k/t)&&i<=x;++i){
                int cnt=0;
                while(x%i==0) ++cnt,x/=i;
                if(cnt>1) goto Next;
                if(cnt) sta[++top]=i,now=-now;
            }
            if(x>1) now=-now,sta[++top]=x;
            for(ri i=0;i<=top;++i) for(ri j=0;j<=Getphi::cnt;++j) f[i][j]=-1;
            if(now== 1) res=add(res,g(n/(k/t),top)*(k/t)%mod);
            if(now==-1) res=dec(res,g(n/(k/t),top)*(k/t)%mod);
            Next:;
        }
    }
    return res*inv%mod;
}
ll n,k1,k2;
int main(){
    // freopen("rand.in","r",stdin);
    // freopen("1.out","w",stdout);
    scanf("%lld%lld%lld",&n,&k1,&k2);
    Getphi::init(n);
    printf("%lld\n",S(n,k1)*S(n,k2)%mod);
}

标签：P7504,frac,R1,int,sum,return,mu,prim,HMOI
来源： https://www.cnblogs.com/Guts/p/15505683.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。