BZOJ 3309: DZY Loves Math

Description

对于正整数$n$，定义$f(n)$为$n$所含质因子的最大幂指数。例如$f(1960)=f(2^3 * 5^1 * 7^2)=3$,$f(10007)=1$,$f(1)=0$
给定正整数$n,m$，求$\sum_{i=1}^n \sum_{j=1}^m f(gcd(i,j))$

Solution

以为自己要独立完成人生中第一道有点难度的数论了。。然而因为常数大被卡掉了。。
化解的过程蛮简单

$$\sum_{i=1}^n \sum_{j=1}^m f(gcd(i,j))$$
$$=\sum_d \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} f(d)*[gcd(i,j)=1]$$
$$=\sum_d \sum_{i=1}^{\lfloor \frac{n}{d}\rfloor} \sum_{j=1}^{\lfloor\frac{m}{d}\rfloor} f(d)*\sum_{d^{'}|gcd(i,j)}\mu(d^{'})$$
令 $a=\lfloor \frac n d \rfloor\ b=\lfloor \frac m d \rfloor$
$$=\sum_d \sum_{i=1}^{a} \sum_{j=1}^{b} f(d)*\sum_{d^{'}}\lfloor\frac a {d^{'}}\rfloor \lfloor \frac b {d^{'}}\rfloor \mu(d^{'})$$
$$=\sum_D\lfloor\frac n D \rfloor \lfloor \frac m D\rfloor \sum_{d|D} f(\frac D d)*\mu(d)$$
然后令 $g(n)=\sum_{d|n}\mu(d)*f(\frac n d)$
原式可以化简为 $$=\sum_D\lfloor\frac n D \rfloor \lfloor \frac m D\rfloor g(D)$$
只要得到 $g$的值，然后按值域封装即可
关键是$g$的值如何得到

我一开始是这么想的：
这种类似卷积形式的求和式都可以以$n ln n$的复杂度得到
我枚举$d$，枚举$\frac n d$，得到$f$
再做一次得到$g$
然后好了
本地测评预处理4S
交上去被卡掉了。。
但蛮有意义放下代码

Code 1

/**************************************************************
    Problem: 3309
    User: bblss123
    Language: C++
    Result: Time_Limit_Exceed
****************************************************************/

#include<iostream>
#include<algorithm>
#include<string.h>
#include<stdio.h>
using namespace std;
typedef long long ll;
const int M=1e7+5;
int _,t,prime[M],g[M],f[M],mu[M];
bool mark[M];
ll sum[M];
inline void Max(int &a,int &b){
    if(a<b)a=b;
}
inline void Euler(){
    mu[1]=1;
    for(int i=2;i<M;++i){
        if(!mark[i])prime[t++]=i,mu[i]=-1;
        for(int j=0,p;j<t&&(p=prime[j])*i<M;++j){
            mark[p*i]=1;
            if(i%p==0){mu[i*p]=0;break;}
            mu[i*p]=-mu[i];
        }
    }
    for(int i=0,d;d=prime[i],i<t;++i){
        for(int k=1,w;(w=k*d)<M;++k)
            Max(f[w],g[w]=g[k]+1);
        for(int k=d;k<M;k+=d)g[k]=0;
    }
    for(int d=1;d<M;++d)
        for(int k=1,w;(w=k*d)<M;++k)
            g[w]+=mu[k]*f[d];
    for(int i=1;i<M;++i)
        sum[i]=sum[i-1]+g[i];
}
inline void rd(int &a){
    a=0;char c;
    while(c=getchar(),!isdigit(c));
    do a=a*10+(c^48);
        while(c=getchar(),isdigit(c));
}
int n,m;
inline void gao(){
    rd(n),rd(m);
    ll ans=0;
    for(int D=1;D<=min(n,m);++D){
        int last=min(n/(n/D),m/(m/D));
        ans+=1ll*(n/D)*(m/D)*(sum[last]-sum[D-1]);
        D=last;
    }
    printf("%lld\n",ans);
}
int main(){
    Euler();
    for(cin>>_;_--;)gao();
}

于是就当我以为又要卡常数的时候，手贱把$g$的值给打出来了，发现值域和$\mu$完全一样。。
然后我开始用容斥推导
意识流地得到了通项式。。
然而意识流，所以去查证明
Claris好神，一眼流地得出通项式，orz 然而等于没看
其他地方连一眼流都没有。。
逛来逛去又回到了PoPoQQQ大爷那
证法通俗易懂，我来抄袭一下

先说公式:
令$n=\Pi_{i=1}^s p_i^{\alpha_i}$
(1):当$\alpha_i$的值完全一样的时候，$g(n)$有非零值$(-1)^{k+1}$
(2):其余时候值为0

证明：
令$mx=\max_{i=1}^s \alpha_i$
注意$\mu(d)$只在$d=\Pi_{i=1}^s p_i^{0\ or\ 1}$时有非零值
那么所有$\alpha_i$至多减$1$
那么我们可以把所有$\alpha_i$划分为两个集合$A,B$
其中$A$中的元素全部等于$mx$，$B$中的元素不等于$mx$
考虑$d$使哪些$\alpha_i$减了$1$
显然$g$值由$A$中的元素决定
假定$A$中有元素被选到
则对于$B$集合的奇偶选择可能性相同，又由$\mu$的性质，知道当前答案为$0$
那么对于$A$中没有元素选取到呢？
同样的，对$B$进行如上讨论，考虑把$B$集合给划分成$A^{'}$和$B^{'}$，然后同样讨论(就好像一个递归)，最后显然可以得到0的答案，则(2)式成立。
现在再来看(1)式就简单多了，直接容斥即可

这样的话 g <script type="math/tex" id="MathJax-Element-344">g</script>可以直接在线性筛里得到

Code 2

/**************************************************************
    Problem: 3309
    User: bblss123
    Language: C++
    Result: Accepted
    Time:11156 ms
    Memory:216136 kb
****************************************************************/

#include<iostream>
#include<algorithm>
#include<string.h>
#include<stdio.h>
using namespace std;
typedef long long ll;
const int M=1e7+5;
int _,t,prime[M>>2],g[M],a[M],pre[M];
ll sum[M];
bool mark[M];
inline void Max(int &a,int &b){
    if(a<b)a=b;
}
inline void Euler(){
    for(int i=2;i<M;++i){
        if(!mark[i])prime[t++]=pre[i]=i,g[i]=a[i]=1;
        for(int j=0,p;j<t&&(p=prime[j])*i<M;++j){
            mark[p*i]=1;
            if(i%p==0){
                a[p*i]=a[i]+1;
                pre[p*i]=pre[i]*p;
                int k=i/pre[i];
                if(k==1)g[p*i]=1;
                else g[p*i]=a[k]==a[p*i]?-g[k]:0;
                break;
            }
            a[i*p]=1;
            pre[i*p]=p;
            if(a[i]==1)g[i*p]=-g[i];
        }
    }
    for(int i=1;i<M;++i)
        sum[i]=sum[i-1]+g[i];
}
inline void rd(int &a){
    a=0;char c;
    while(c=getchar(),!isdigit(c));
    do a=a*10+(c^48);
        while(c=getchar(),isdigit(c));
}
int n,m;
inline int gcd(int n,int m){
    return m?gcd(m,n%m):n;
}
inline void nt(ll x){
    if(!x)return;
    nt(x/10);
    putchar(48+x%10);
}
inline void pt(ll x){
    if(!x)putchar('0');
    else nt(x);
}
inline void gao(){
    rd(n),rd(m);
    ll ans=0;
    for(int D=1;D<=min(n,m);++D){
        int last=min(n/(n/D),m/(m/D));
        ans+=1ll*(n/D)*(m/D)*(sum[last]-sum[D-1]);
        D=last;
    }
    pt(ans),putchar('\n');
}
int main(){
    Euler();
    for(rd(_);_--;)gao();
}