P7504 「HMOI R1」可愛的德麗莎

...

好像被我推複雜了...

簡述題意

求：

\[\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]\)。

那麼有：

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))\\

如果是求 \(i\varphi(i)\) 的字首和，可以直接上

Min_25

，可是這裡的形式比較詭異，不能直接求。

不過注意到 \(\mu(t)\not = 0\)，當且僅當 \(t\) 中的每一個質因子隻出現一次，似乎是可以做一個DP。

我們令 \(p\) 表示 \(t\) 中的某一個質因子，同時定義一個新的函數 \(G(n,t)\)。

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)

加上記憶化之後，複雜度便正确了。

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