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.csdn.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 }
upd181010:
$sum_{k|T}kmu(frac{T}{k})(frac{T}{k})^2=Tsum_{k|T}mu(frac{T}{k})frac{T}{k}=Tsum_{k|T}mu(k)k$
这个式子似乎并没有什么用...但是可以记一下，说不定下次就给出最后的那个呢
bzoj3309
https://www.lydsy.com/JudgeOnline/problem.php?id=3309
不妨设a<=b
显然原式=$sum_{k=1}^bf(k)sum_{i=1}^{{lfloor}frac{b}{k}{ floor}}mu(i){{lfloor}frac{a}{ik}{ floor}}{{lfloor}frac{b}{ik}{ floor}}$
用跟上面类似的方法变到$sum_{T=1}^b{lfloor}frac{a}{T}{ floor}{lfloor}frac{b}{T}{ floor}sum_{k|T}f(k)mu(frac{T}{k})$
后面那个不太好筛啊，根本就不会，我又去查题解了。。。
（弃用）如果i%p==0，可以丢掉那些新产生因子的$mu$，每一项是由$mu(d)f(i/d)$变到$mu(d)f(i/d*p)$，显然会加(i的(i/d)中"最大幂指数等于最小质因子的幂指数"的数的$mu(d)$之和)，那个求时依赖的值也可以递推出
（弃用）如果i%p!=0，那么首先会多一些产生贡献的项，这些项是由原来i的任意因子乘上一个p得到的，是$mu(i/d)f(d)$变到$mu(i*p/(d*p))f(d*p)=mu(i/d)f(dp)$，产生的贡献是原来的总和加上一个值（i的因子中"最大幂指数等于最小质因子的幂指数"的数(设为d)的$mu(i/d)$之和）；
筛法的解释：https://blog.csdn.net/phenix_2015/article/details/50799021
这条结论，简单来讲，就是那个函数，只有各个质因子的指数相等的有值((-1)^{质因子数+1})，其他都为0
直接想好像很难想到的样子。。以后记得做这类题要打打表分解分解质因数找找规律
算是A了吧。。。
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 N 10000000 15 ll prime[N+100],len,mu[N+100]; 16 ll n1[N+100],n2[N+100],n3[N+100],h[N+100]; 17 //分别为最小质因子次数，其他质因子次数，质因子数，函数值 18 //n2为-1:第一个质因子;n2为-2:已经失败 19 bool nprime[N+100]; 20 ll a,b,ans,T; 21 template<class T> 22 inline void read(T &x) { 23 int f=1;x=0;char ch=getchar(); 24 while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();} 25 while(ch>='0'&&ch<='9'){x*=10;x+=(ch-'0');ch=getchar();} 26 x*=f; 27 } 28 template<class T> 29 inline void write(T x) { 30 if(x<0) putchar('-'),x=-x; 31 if(x>9) write(x/10); 32 putchar(x%10+'0'); 33 } 34 int main() 35 { 36 ll i,j; 37 mu[1]=1; 38 for(i=2;i<=N;i++) 39 { 40 if(!nprime[i]) 41 { 42 prime[++len]=i; 43 mu[i]=-1; 44 n2[i]=-1; 45 n3[i]=1; 46 n1[i]=1; 47 h[i]=1; 48 } 49 for(j=1;j<=len&&i*prime[j]<=N;j++) 50 { 51 nprime[i*prime[j]]=1; 52 if(i%prime[j]==0) 53 { 54 mu[i*prime[j]]=0; 55 n1[i*prime[j]]=n1[i]+1; 56 n2[i*prime[j]]=n2[i]; 57 n3[i*prime[j]]=n3[i]; 58 if(n2[i*prime[j]]==-1||n1[i*prime[j]]==n2[i*prime[j]]) 59 h[i*prime[j]]=(n3[i*prime[j]]%2==0)?-1:1; 60 else 61 h[i*prime[j]]=0; 62 break; 63 } 64 else 65 { 66 mu[i*prime[j]]=-mu[i]; 67 n1[i*prime[j]]=1; 68 n2[i*prime[j]]=(n2[i]==-1||n2[i]==n1[i])?n1[i]:-2; 69 n3[i*prime[j]]=n3[i]+1; 70 if(n2[i*prime[j]]==-1||n1[i*prime[j]]==n2[i*prime[j]]) 71 h[i*prime[j]]=(n3[i*prime[j]]%2==0)?-1:1; 72 else 73 h[i*prime[j]]=0; 74 } 75 } 76 } 77 //for(i=1;i<=1000;i++) printf("%lld %lld %lld %lld ",i,h[i],n1[i],n2[i]); 78 //return 0; 79 for(i=1;i<=N;i++) h[i]+=h[i-1]; 80 read(T); 81 //T=1000; 82 while(T--) 83 { 84 read(a);read(b); 85 //a=7558588;b=9653114; 86 if(a==0||b==0) 87 { 88 puts("0"); 89 continue; 90 } 91 if(a>b) swap(a,b); 92 ans=0; 93 for(i=1;i<=a;i=j+1) 94 { 95 j=min(a/(a/i),b/(b/i)); 96 ans+=(a/i)*(b/i)*(h[j]-h[i-1]); 97 } 98 write(ans);putchar(' '); 99 } 100 return 0; 101 }