bzoj2154||洛谷P1829 Crash的数字表格&&JZPTAB && bzoj3309 DZY Loves Math
bzoj2154||洛谷P1829
https://www.lydsy.com/JudgeOnline/problem.php?id=2154
https://www.luogu.org/problemnew/show/P1829
不妨设n<=m
就是求$ans=sum_{k=1}^m{frac{1}{k}}sum_{i=1}^nsum_{j=1}^m{ij[(i,j)=k]}$
把1/k后面的那一部分提出来，设为f(k)，
然后莫比乌斯反演得到f(k)较简易的计算式，代回ans，并化简（过程略）
最后化简出来是$sum_{k=1}^m{k}sum_{i=1}^{{lfloor}frac{m}{k}{ floor}}mu(i){i^2}g({{lfloor}frac{{lfloor}frac{n}{k}{ floor}}{i}{ floor}},{{lfloor}frac{{lfloor}frac{m}{k}{ floor}}{i}{ floor}})$
其中$g(x,y)=frac{(x+1)x(y+1)y}{4}$
可以通过两重的整除分块以O(n)的复杂度计算（第一重枚举n/k,m/k，第二重枚举(n/k)/i,(m/k)/i）
1 #pragma GCC optimize(3) 2 #include<cstdio> 3 #include<algorithm> 4 #include<cstring> 5 #include<vector> 6 using namespace std; 7 #define fi first 8 #define se second 9 #define mp make_pair 10 #define pb push_back 11 typedef long long ll; 12 typedef unsigned long long ull; 13 typedef pair<int,int> pii; 14 #define md 20101009 15 #define N 10000000 16 ll prime[N+100],len,mu[N+100],dd[N+100]; 17 bool nprime[N+100]; 18 ll Mod(ll a,ll b) 19 { 20 if(a>=0) return a%b; 21 else if(a%b==0) return 0; 22 else return b+a%b; 23 } 24 ll G(ll x,ll y) 25 { 26 return (x+1)*x%md*(y+1)%md*y%md*15075757%md; 27 } 28 ll calc(ll n,ll m) 29 { 30 ll i,j,ans=0; 31 for(i=1;i<=m;i=j+1) 32 { 33 if(n<i) j=min(m,m/(m/i)); 34 else j=min(m,min(n/(n/i),m/(m/i))); 35 ans=Mod(ans+(dd[j]-dd[i-1])*G(n/i,m/i),md); 36 } 37 return ans; 38 } 39 ll n,m; 40 int main() 41 { 42 ll i,j,ans=0; 43 mu[1]=1; 44 for(i=2;i<=N;i++) 45 { 46 if(!nprime[i]) prime[++len]=i,mu[i]=-1; 47 for(j=1;j<=len&&i*prime[j]<=N;j++) 48 { 49 nprime[i*prime[j]]=1; 50 if(i%prime[j]==0) {mu[i*prime[j]]=0;break;} 51 else mu[i*prime[j]]=-mu[i]; 52 } 53 } 54 for(i=1;i<=N;i++) dd[i]=dd[i-1]+i*i*mu[i]; 55 scanf("%lld%lld",&n,&m); 56 //n=10000000;m=1000; 57 if(n>m) swap(n,m); 58 for(i=1;i<=m;i=j+1) 59 { 60 if(n<i) j=min(m,m/(m/i)); 61 else j=min(m,min(n/(n/i),m/(m/i))); 62 ans=(ans+(i+j)*(j-i+1)/2*calc(n/i,m/i))%md; 63 } 64 printf("%lld",ans); 65 return 0; 66 }
额，好像还有一道题面基本一样的题（然而bzoj权限题，不能交），然而是很多组数据，不能每次O(n)算
相关：https://blog.****.net/qq_30974369/article/details/79087445
对于这个式子$sum_{k=1}^m{k}sum_{i=1}^{{lfloor}frac{m}{k}{ floor}}mu(i){i^2}g({{lfloor}frac{n}{ik}{ floor}},{{lfloor}frac{m}{ik}{ floor}})$
令T=ik
原式=$sum_{k=1}^m{k}sum_{k|T}mu(frac{T}{k})(frac{T}{k})^2g({{lfloor}frac{n}{T}{ floor}},{{lfloor}frac{m}{T}{ floor}})$
=$sum_{T=1}^mg({{lfloor}frac{n}{T}{ floor}},{{lfloor}frac{m}{T}{ floor}})sum_{k|T}kmu(frac{T}{k})(frac{T}{k})^2$
最后从sigma k|T开始那一部分，等于$(id*(mucdot idcdot id))(T)$
是个积性函数，且是可以线性筛出来的（然而我不会，又去查了题解），将其设为函数h
对于1，h[1]=1
对于一个质数p，h[p]=$1*mu(p)*p^2+p*mu(1)*1^2$
对于i%p!=0，h(i*p)=h(i)*h(p)（由于积性）
对于i%p==0，由于h(i*p)拆开后式子中比h(i)多的项的$mu$值均为0，因此只有对于h(i)中已有的项增加贡献；每一项由$mu(d)*d^2*(i/d)$变到$mu(d)*d^2*((i*p)/d)$，因此h(i*p)=h(i)*p
筛出来之后做一下前缀和，这样子每一次询问就可以数论分块根号解决了
1 #include<cstdio> 2 #include<algorithm> 3 #include<cstring> 4 #include<vector> 5 using namespace std; 6 #define fi first 7 #define se second 8 #define mp make_pair 9 #define pb push_back 10 typedef long long ll; 11 typedef unsigned long long ull; 12 typedef pair<int,int> pii; 13 #define md 20101009 14 #define N 10000000 15 ll prime[N+100],len,dd[N+100]; 16 bool nprime[N+100]; 17 ll h[N+100]; 18 ll Mod(ll a,ll b) 19 { 20 if(a>=0) return a%b; 21 else if(a%b==0) return 0; 22 else return b+a%b; 23 } 24 ll G(ll x,ll y) 25 { 26 return (x+1)*x%md*(y+1)%md*y%md*15075757%md; 27 } 28 ll n,m; 29 int main() 30 { 31 ll i,j,ans=0; 32 h[1]=1; 33 for(i=2;i<=N;i++) 34 { 35 if(!nprime[i]) prime[++len]=i,h[i]=Mod(-i*i+i,md); 36 for(j=1;j<=len&&i*prime[j]<=N;j++) 37 { 38 nprime[i*prime[j]]=1; 39 if(i%prime[j]==0) {h[i*prime[j]]=h[i]*prime[j]%md;break;} 40 else h[i*prime[j]]=h[i]*h[prime[j]]%md; 41 } 42 } 43 for(i=1;i<=N;i++) dd[i]=(dd[i-1]+h[i])%md; 44 scanf("%lld%lld",&n,&m); 45 //n=10000000;m=1000; 46 if(n>m) swap(n,m); 47 for(i=1;i<=m;i=j+1) 48 { 49 if(n<i) j=min(m,m/(m/i)); 50 else j=min(m,min(n/(n/i),m/(m/i))); 51 ans=(ans+Mod(dd[j]-dd[i-1],md)*G(n/i,m/i))%md; 52 } 53 printf("%lld",ans); 54 return 0; 55 }