![[国家集训队]Crash的数字表格/JZPTAB和51Nod1238最小公倍数之和V3题解](/uploads/image/0119.jpg)
给你两个正整数 n , m n,m n,m，求下面式子在 m o d    20101009 mod 20101009 mod20101009意义下的值。
∑ i = 1 n ∑ j = 1 m l c m ( i , j ) sum_{i=1}^nsum_{j=1}^mlcm(i,j) i=1∑nj=1∑mlcm(i,j)
n , m ≤ 1 0 7 n,mleq 10^7 n,m≤107
我们根据 l c m ( i j ) = i j g c d ( i , j ) lcm(ij)=frac{ij}{gcd(i,j)} lcm(ij)=gcd(i,j)ij，将原式化简成：
∑ i = 1 n ∑ j = 1 m i j g c d ( i , j ) sum_{i=1}^nsum_{j=1}^mfrac{ij}{gcd(i,j)} i=1∑nj=1∑mgcd(i,j)ij
根据套路，我们枚举 g c d gcd gcd，原式则化成：
（这里默认 n ≤ m nleq m n≤m，如果不满足则交换）
= ∑ d = 1 n ∑ i = 1 n ∑ j = 1 m i j d [ g c d ( i , j ) = d ] = ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ ( i d ) × ( j d ) d [ g c d ( i , j ) = 1 ] =sum_{d=1}^nsum_{i=1}^nsum_{j=1}^mfrac{ij}{d}[gcd(i,j)=d]\ = sum_{d=1}^nsum_{i=1}^{lfloorfrac{n}{d}rfloor}sum_{j=1}^{lfloorfrac{m}{d}rfloor}frac{(id)times(jd)}{d}[gcd(i,j)=1] =d=1∑ni=1∑nj=1∑mdij[gcd(i,j)=d]=d=1∑ni=1∑⌊dn⌋j=1∑⌊dm⌋d(id)×(jd)[gcd(i,j)=1]
此时原式就可以变成：
∑ d = 1 n d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ i j [ g c d ( i , j ) = 1 ] sum_{d=1}^ndsum_{i=1}^{lfloorfrac{n}{d}rfloor}sum_{j=1}^{lfloorfrac{m}{d}rfloor}ij[gcd(i,j)=1] d=1∑ndi=1∑⌊dn⌋j=1∑⌊dm⌋ij[gcd(i,j)=1]
然后又用莫比乌斯反演，我们可以将原式变成：
= ∑ d = 1 n d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ m d ⌋ i j ∑ w ∣ i , w ∣ j μ ( w ) = ∑ d = 1 n d ∑ w = 1 ⌊ n d ⌋ μ ( w ) ∑ i = 1 ⌊ n d w ⌋ i w ∑ j = 1 ⌊ n d w ⌋ j w =sum_{d=1}^ndsum_{i=1}^{lfloorfrac{n}{d}rfloor}sum_{j=1}^{lfloorfrac{m}{d}rfloor}ijsum_{w|i,w|j}mu(w)\ =sum_{d=1}^ndsum_{w=1}^{lfloorfrac{n}{d}rfloor}mu(w)sum_{i=1}^{lfloorfrac{n}{dw}rfloor}iwsum_{j=1}^{lfloorfrac{n}{dw}rfloor}jw =d=1∑ndi=1∑⌊dn⌋j=1∑⌊dm⌋ijw∣i,w∣j∑μ(w)=d=1∑ndw=1∑⌊dn⌋μ(w)i=1∑⌊dwn⌋iwj=1∑⌊dwn⌋jw
然后我们令 S ( n ) = ∑ i = 1 n i = n × ( n + 1 ) 2 S(n)=sum_{i=1}^ni=frac{ntimes(n+1)}{2} S(n)=∑i=1ni=2n×(n+1)，原式就可以变成：
= ∑ d = 1 n d ∑ w = 1 ⌊ n d ⌋ μ ( w ) w 2 S ( ⌊ n d w ⌋ ) S ( ⌊ m d w ⌋ ) =sum_{d=1}^ndsum_{w=1}^{lfloorfrac{n}{d}rfloor}mu(w)w^2S(leftlfloorfrac{n}{dw}rightrfloor)S(leftlfloorfrac{m}{dw}rightrfloor) =d=1∑ndw=1∑⌊dn⌋μ(w)w2S(⌊dwn⌋)S(⌊dwm⌋)
用线性筛筛出 μ ( w ) w 2 mu(w)w^2 μ(w)w2的前缀和，然后 S S S可以 O ( 1 ) O(1) O(1)算，所以我们可以对于 H ( n , m ) = ∑ w = 1 n μ ( w ) w 2 S ( ⌊ n w ⌋ ) S ( ⌊ m w ⌋ ) H(n,m)=sum_{w=1}^{n}mu(w)w^2S(leftlfloorfrac{n}{w}rightrfloor)S(leftlfloorfrac{m}{w}rightrfloor) H(n,m)=∑w=1nμ(w)w2S(⌊wn⌋)S(⌊wm⌋)用数论分块计算，而对于 ∑ d = 1 n H ( ⌊ n d ⌋ , ⌊ m d ⌋ ) sum_{d=1}^nH(leftlfloorfrac{n}{d}rightrfloor,leftlfloorfrac{m}{d}rightrfloor) ∑d=1nH(⌊dn⌋,⌊dm⌋)又数论分块算即可，复杂度最坏为 O ( n + m ) O(n+m) O(n+m)。
还是上面的问题， n , m ≤ 1 0 7 n,mleq10^7 n,m≤107，但是有 T ≤ 1 0 4 Tleq 10^4 T≤104组询问。
所以我们必须要做到 O ( n ) O(n) O(n)预处理，每次 O ( n ) O(sqrt{n}) O(n )的回答。
我们还是看上一个化简后的式子：
= ∑ d = 1 n d ∑ w = 1 ⌊ n d ⌋ μ ( w ) w 2 S ( ⌊ n d w ⌋ ) S ( ⌊ m d w ⌋ ) =sum_{d=1}^ndsum_{w=1}^{lfloorfrac{n}{d}rfloor}mu(w)w^2S(leftlfloorfrac{n}{dw}rightrfloor)S(leftlfloorfrac{m}{dw}rightrfloor) =d=1∑ndw=1∑⌊dn⌋μ(w)w2S(⌊dwn⌋)S(⌊dwm⌋)
还有一个套路就是，我们枚举 d w dw dw，所以令 T = d w T=dw T=dw，那么原式等于：
= ∑ T = 1 n S ( ⌊ n T ⌋ ) S ( ⌊ m T ⌋ ) ∑ d ∣ T d μ ( T d ) ( T d ) 2 = ∑ T = 1 n S ( ⌊ n T ⌋ ) S ( ⌊ m T ⌋ ) ∑ w ∣ T μ ( w ) w 2 T w = ∑ T = 1 n S ( ⌊ n T ⌋ ) S ( ⌊ m T ⌋ ) T ∑ w ∣ T μ ( w ) w =sum_{T=1}^nS(leftlfloorfrac{n}{T}rightrfloor)S(leftlfloorfrac{m}{T}rightrfloor)sum_{d|T}dmu(frac{T}{d})(frac{T}{d})^2\ =sum_{T=1}^nS(leftlfloorfrac{n}{T}rightrfloor)S(leftlfloorfrac{m}{T}rightrfloor)sum_{w|T}mu(w)w^2frac{T}{w}\ =sum_{T=1}^nS(leftlfloorfrac{n}{T}rightrfloor)S(leftlfloorfrac{m}{T}rightrfloor)Tsum_{w|T}mu(w)w =T=1∑nS(⌊Tn⌋)S(⌊Tm⌋)d∣T∑dμ(dT)(dT)2=T=1∑nS(⌊Tn⌋)S(⌊Tm⌋)w∣T∑μ(w)w2wT=T=1∑nS(⌊Tn⌋)S(⌊Tm⌋)Tw∣T∑μ(w)w
对于前半部分我们每次数论分块求即可，而对于后面我们线性筛出它，然后求前缀和即可。
我们令 f ( x ) = x ∑ d ∣ x μ ( x ) x f(x)=xsum_{d|x}mu(x)x f(x)=x∑d∣xμ(x)x，那么它肯定是个积性函数（因为相当于可以写成 i d ⋅ ( ( μ ⋅ i d ) ⨂ 1 ) idcdot ((mucdot id)bigotimes 1) id⋅((μ⋅id)⨂1)）根据欧拉筛的方法我们只需维护三个东西 c n t , f i r , p o w cnt,fir,pow cnt,fir,pow即可，分别表示： c n t [ i ] cnt[i] cnt[i]， i i i的最小质因子指数； f i r [ i ] fir[i] fir[i]， i i i的最小质因子； p o w [ i ] pow[i] pow[i]表示 f i r [ i ] c n t [ i ] fir[i]^{cnt[i]} fir[i]cnt[i]。
然后我们来看 f ( x ) f(x) f(x)的计算方法：
代码就十分好写了，同样也可以过第一题：
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int M=1e7+10,P=1e6+10;
const ll Mod=20101009;
ll prime[P],cnt;bool vis[M];
ll F[M],c[M],p[M],f[M],inv_2;
ll fpow(ll a,ll b){ll ans=1;for(;b;b>>=1,a=(a*a)%Mod){if(b&1)ans=(ans*a)%Mod;}return ans;
}
void init(ll n){F[1]=1;for(int i=1;i<=n;i++)p[i]=f[i]=1;for(ll i=2;i<=n;i++){if(!vis[i]){prime[++cnt]=i;F[i]=(i-(i*i%Mod)+Mod)%Mod;c[i]=1;f[i]=i;p[i]=i;}for(ll j=1,v;j<=cnt&&i*prime[j]<=n;j++){v=i*prime[j];vis[v]=1;if(!(i%prime[j])){c[v]=c[i]+1;f[v]=f[i];p[v]=p[i]*f[i];F[v]=F[v/p[v]]*((p[v]-p[v]*f[v]%Mod+Mod)%Mod)%Mod;break;}F[v]=(F[i]*F[prime[j]])%Mod;c[v]=1;f[v]=prime[j];p[v]=prime[j];}}for(int i=2;i<=n;i++)F[i]=(F[i]+F[i-1])%Mod;
}
ll S(ll a){if(a>=Mod)a%=Mod;return ((a*(a+1)%Mod)*inv_2)%Mod;}
ll calc(ll L,ll R){return ((F[R]-F[L-1])%Mod+Mod)%Mod;}
ll solve(ll n,ll m){if(n>m)swap(n,m);init(n);inv_2=fpow(2,Mod-2);ll ans=0;for(ll i=1,j;i<=n;i=j+1){j=min(n/(n/i),m/(m/i));ans=(ans+S(n/i)*S(m/i)%Mod*calc(i,j)%Mod)%Mod;}return ans;
}
ll n,m;
int main(){scanf("%lld%lld",&n,&m);printf("%lldn",solve(n,m));return 0;
}
还是一样的问题，只不过 n = m ≤ 1 0 10 n=mleq10^{10} n=m≤1010，只询问一次。
此时就不能线性筛了，必须要做到线性以下的求法，那么就是杜教筛了。
而前面化简的式子，若 f ( x ) = x ∑ d ∣ x μ ( d ) d f(x)=xsum_{d|x}mu(d)d f(x)=x∑d∣xμ(d)d能杜教筛出，那么就可以用之前的方法做，但是这并不好实现反正博主不会的啦QWQ。
所以我们考虑，在这个时候重新变换：
= ∑ d = 1 n d ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 ⌊ n d ⌋ i j [ g c d ( i , j ) = 1 ] =sum_{d=1}^ndsum_{i=1}^{lfloorfrac{n}{d}rfloor}sum_{j=1}^{lfloorfrac{n}{d}rfloor}ij[gcd(i,j)=1] =d=1∑ndi=1∑⌊dn⌋j=1∑⌊dn⌋ij[gcd(i,j)=1]
我们可以先来看个式子：
∑ i = 1 n ∑ j = 1 i [ g c d ( i , j ) = 1 ] i j sum_{i=1}^{n}sum_{j=1}^i[gcd(i,j)=1]ij i=1∑nj=1∑i[gcd(i,j)=1]ij
我们知道 φ ( i ) varphi(i) φ(i)表示 1 ∼ i − 1 1sim i-1 1∼i−1中与 i i i互质的个数，而容易得知， 1 ∼ i − 1 1sim i-1 1∼i−1中与 i − 1 i-1 i−1互质的数的和为 i φ ( i ) 2 frac{ivarphi(i)}{2} 2iφ(i)
简单证明:因为对于一个数 d ≤ i − 1 dleq i-1 d≤i−1，如果它与 i i i互质，根据欧几里得算法， i − d i-d i−d也与 i i i互质，那么这一对互质的和就为 i i i，总共有 φ ( i ) varphi(i) φ(i)个互质，总和就是 i φ ( i ) ivarphi(i) iφ(i)，但是 d , i − d d,i-d d,i−d和 i − d , d i-d,d i−d,d算了两次，所以再除以2就是答案了，其中1要特殊考虑。
所以上面的那个式子可以写成：
= ∑ i = 1 n i ∑ j = 1 i j [ g c d ( i , j ) = 1 ] = ∑ i = 1 n i i φ ( i ) + [ i = 1 ] 2 =sum_{i=1}^nisum_{j=1}^ij[gcd(i,j)=1]\ =sum_{i=1}^nifrac{ivarphi(i)+[i=1]}{2} =i=1∑nij=1∑ij[gcd(i,j)=1]=i=1∑ni2iφ(i)+[i=1]
那么我们看原式，就可以写成：
= ∑ d = 1 n d ( 2 × ( ∑ i = 1 ⌊ n d ⌋ i i φ ( i ) + [ i = 1 ] 2 ) − 1 ) = ∑ d = 1 n d ( ∑ i = 1 ⌊ n d ⌋ i 2 φ ( i ) ) =sum_{d=1}^ndleft(2timesleft(sum_{i=1}^{lfloorfrac{n}{d}rfloor}ifrac{ivarphi(i)+[i=1]}{2}right)-1right)\ =sum_{d=1}^ndleft(sum_{i=1}^{lfloorfrac{n}{d}rfloor}{i^2varphi(i)}right) =d=1∑nd⎝⎛2×⎝⎛i=1∑⌊dn⌋i2iφ(i)+[i=1]⎠⎞−1⎠⎞=d=1∑nd⎝⎛i=1∑⌊dn⌋i2φ(i)⎠⎞
因为原式后半部分可以写成：
∑ i = 1 ⌊ n d ⌋ i ( ∑ j = 1 i j [ g c d ( i , j ) = 1 ] + ∑ j = i + 1 ⌊ n d ⌋ j [ g c d ( i , j ) = 1 ] ) sum_{i=1}^{lfloorfrac{n}{d}rfloor}ileft(sum_{j=1}^{i}j[gcd(i,j)=1]+sum_{j=i+1}^{lfloorfrac{n}{d}rfloor}j[gcd(i,j)=1]right) i=1∑⌊dn⌋i⎝⎛j=1∑ij[gcd(i,j)=1]+j=i+1∑⌊dn⌋j[gcd(i,j)=1]⎠⎞
而后面两部分就想一个正方形，边长为 ⌊ n d ⌋ lfloorfrac{n}{d}rfloor ⌊dn⌋，沿着对角线剪开一样，是对称的，所以取其中一个乘以 2 2 2即可。
详细来说，原式可以变成：
∑ d = 1 n d ( 2 × ( ∑ i = 1 ⌊ n d ⌋ ∑ j = 1 i i j ) − ∑ i = 1 ⌊ n d ⌋ i 2 ) sum_{d=1}^nd left(2times left(sum_{i=1}^{lfloorfrac{n}{d}rfloor}sum_{j=1}^iijright)-sum_{i=1}^{lfloorfrac{n}{d}rfloor}i^2right) d=1∑nd⎝⎛2×⎝⎛i=1∑⌊dn⌋j=1∑iij⎠⎞−i=1∑⌊dn⌋i2⎠⎞
因为 ( i , j ) , ( j , i ) (i,j),(j,i) (i,j),(j,i)会被算两次，所以乘以2，而 ( i , i ) (i,i) (i,i)只会被算一次所以减去，然后就可以用前面式子替换就可以得到，后面 ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ i 2 = n × ( n + 1 ) 2 sum_{d=1}^nsum_{i=1}^{lfloorfrac{n}{d}rfloor}i^2=frac{ntimes(n+1)}{2} ∑d=1n∑i=1⌊dn⌋i2=2n×(n+1)减去 ∑ d = 1 n d ∑ i = 1 ⌊ n d ⌋ [ i = 1 ] = n × ( n + 1 ) 2 sum_{d=1}^ndsum_{i=1}^{lfloorfrac{n}{d}rfloor}[i=1]=frac{ntimes(n+1)}{2} ∑d=1nd∑i=1⌊dn⌋[i=1]=2n×(n+1)的刚好消掉。
那么原式就转换成了求：
∑ d = 1 n d ( ∑ i = 1 ⌊ n d ⌋ i 2 φ ( i ) ) sum_{d=1}^ndleft(sum_{i=1}^{lfloorfrac{n}{d}rfloor}{i^2varphi(i)}right) d=1∑nd⎝⎛i=1∑⌊dn⌋i2φ(i)⎠⎞
我们令 f ( x ) = x 2 φ ( x ) f(x)=x^2varphi(x) f(x)=x2φ(x)，然后令 g ( x ) = x 2 g(x)=x^2 g(x)=x2，容易得知（ ⨂ bigotimes ⨂为狄利克雷卷积）：
f ⨂ g ( n ) = ∑ d ∣ n f ( d ) g ( n d ) = ∑ d ∣ n d 2 φ ( d ) ( n d ) 2 = n 2 ∑ d ∣ n φ ( d ) fbigotimes g(n)=sum_{d|n}f(d)g(frac{n}{d})=sum_{d|n}d^2varphi(d)(frac{n}{d})^2=n^2sum_{d|n}varphi(d) f⨂g(n)=d∣n∑f(d)g(dn)=d∣n∑d2φ(d)(dn)2=n2d∣n∑φ(d)
又因为 ∑ d ∣ n φ ( d ) = n sum_{d|n}varphi(d)=n ∑d∣nφ(d)=n，所以令：
h ( n ) = f ⨂ g ( n ) = n 3 h(n)=fbigotimes g(n)=n^3 h(n)=f⨂g(n)=n3
令 S ( n ) = ∑ i = 1 n f ( i ) S(n)=sum_{i=1}^nf(i) S(n)=∑i=1nf(i)，因为 g ( 1 ) = 1 g(1)=1 g(1)=1
带入杜教筛公式中可以得到：
S ( n ) = ∑ i = 1 n h ( i ) − ∑ i = 2 n g ( i ) S ( ⌊ n i ⌋ ) S(n)=sum_{i=1}^nh(i)-sum_{i=2}^ng(i)S(leftlfloorfrac{n}{i}rightrfloor) S(n)=i=1∑nh(i)−i=2∑ng(i)S(⌊in⌋)
也就是：
S ( n ) = ( n × ( n + 1 ) 2 ) 2 − ∑ i = 2 n i 2 S ( ⌊ n i ⌋ ) S(n)=left(frac{ntimes(n+1)}{2}right)^2-sum_{i=2}^ni^2S(leftlfloorfrac{n}{i}rightrfloor) S(n)=(2n×(n+1))2−i=2∑ni2S(⌊in⌋)
而 ∑ i = 1 n i 2 = n × ( n + 1 ) × ( 2 n + 1 ) 6 sum_{i=1}^ni^2=frac{ntimes(n+1)times(2n+1)}{6} ∑i=1ni2=6n×(n+1)×(2n+1)，所以这个式子就很容易算了。
用杜教筛的方式先预处理 n 2 3 n^{frac{2}{3}} n32的答案，然后数论分块递归即可。
对于原式外面我们仍然数论分块求即可。
复杂度 O ( n 2 3 ) O(n^{frac{2}{3}}) O(n32)
#include<cstdio>
#include<cstring>
#include<algorithm>
#define ll long long
using namespace std;
const int N=1e6+10,M=1e4+10;
const ll Mod=1e9+7;
ll prime[N],phi[N],cnt,n;
ll F[N],D[M],inv_2,inv_6;
ll fpow(ll a,ll b){ll ans=1;for(;b;b>>=1,a=(a*a)%Mod){if(b&1)ans=(ans*a)%Mod;}return ans;
}
ll &Rec(ll a){if(a<N) return F[a];else return D[n/a];
}
bool vis[N];
ll Sqr(ll a){if(a>=Mod)a%=Mod;return a*a%Mod;}
void init(int up){memset(F,-1,sizeof(F));memset(D,-1,sizeof(D));phi[1]=1;inv_2=fpow(2,Mod-2);inv_6=fpow(6,Mod-2);for(int i=2;i<up;i++){if(!vis[i]){prime[++cnt]=i;phi[i]=i-1;}for(int j=1,v;j<=cnt&&i*prime[j]<up;j++){v=i*prime[j];vis[v]=1;if(!(i%prime[j])){phi[v]=phi[i]*prime[j];break;}phi[v]=phi[i]*phi[prime[j]];}}F[0]=0;for(int i=1;i<up;i++)F[i]=(F[i-1]+phi[i]*Sqr(i)%Mod)%Mod;
}
ll S(ll a){if(a>=Mod)a%=Mod;return (a*(a+1)%Mod)*inv_2%Mod;}
ll S2(ll a){if(a>=Mod)a%=Mod;return (a*(a+1)%Mod*((a+a+1)%Mod)%Mod)*inv_6%Mod;}
ll S3(ll a){if(a>=Mod)a%=Mod;return Sqr(S(a));}
ll Sarea(ll L,ll R){return ((S(R)-S(L-1))%Mod+Mod)%Mod;}
ll S2area(ll L,ll R){return ((S2(R)-S2(L-1))%Mod+Mod)%Mod;}
ll calc(ll x){ll &ans=Rec(x);if(ans!=-1) return ans;ans=S3(x);for(ll i=2,j;i<=x;i=j+1){j=x/(x/i);(ans-=(S2area(i,j)*calc(x/i))%Mod)%=Mod;}return (ans+Mod)%Mod;
}
ll solve(){ll ans=0;for(ll i=1,j;i<=n;i=j+1){j=n/(n/i);(ans+=(Sarea(i,j)*calc(n/i))%Mod)%=Mod;}return ans;
}
int main(){scanf("%lld",&n);init(min(n+1,1ll*N));printf("%lldn",solve());return 0;
}
大佬的更好的讲解-by fwat Orz%%%【链接】
另一种方法-by Imagine Orz%%%【链接】
如果有错，或者有更好的方法，欢迎提出！(*｀∀´*)ノ亻!
快要退役了，再努力搏一搏吧
本文发布于:2024-02-05 03:31:02，感谢您对本站的认可！
本文链接：https://www.4u4v.net/it/170723039762680.html
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系，我们将在24小时内删除。
| 留言与评论（共有 0 条评论） |
