BZOJ 3512: DZY Loves Math IV [杜教筛]

BZOJ 3512: DZY Loves Math IV [杜教筛]

注意是去除所有质因子的乘积的n

厉害的思想是构造互质情况把phi(i*j)分开

补充candy？的最后一步推导：枚举e，那么gcd(i,n)要是e的倍数，枚举是e的多少倍（上界[m/e]），乘上phi(i*d)，就是S(n,m)的形式了

#include
#define reg register int
#define il inline
#define numb (ch^'0')
using namespace std;
typedef long long ll;
il void rd(int &x){char ch;x=0;bool fl=false;while(!isdigit(ch=getchar()))(ch=='-')&&(fl=true);for(x=numb;isdigit(ch=getchar());x=x*10+numb);(fl==true)&&(x=-x);
}
namespace Miracle{
const int N=1e5+5;
const int M=1999880;
const int mod=1e9+7;
int pri[M],tot;
int pro[M];
int phi[M],vis[M];
int sum[M];
int n,m;
int ad(int x,int y){return (x+y)>=mod?x+y-mod:x+y;
}
void sieve(int n){phi[1]=sum[1]=1;pro[1]=1;for(reg i=2;i<=n;++i){if(!vis[i]){pri[++tot]=i;phi[i]=i-1;pro[i]=i;}for(reg j=1;j<=tot;++j){if(i*pri[j]>n) break;vis[i*pri[j]]=1;if(i%pri[j]==0){phi[i*pri[j]]=phi[i]*pri[j];pro[i*pri[j]]=pro[i];break;}phi[i*pri[j]]=phi[i]*phi[pri[j]];pro[i*pri[j]]=pro[i]*pri[j];}sum[i]=ad(sum[i-1],phi[i]);}
}
const int G=31640;
int p1[G],p2[G];
int Phi(int n){if(n<=M-5){return sum[n];}//cout<<" 23333 ------"<if(n<G){if(p1[n]!=-1) return p1[n];}else{if(p2[m/n]!=-1) return p2[m/n];}int ret=(ll)n*(n+1)/2%mod;for(reg i=2,x=0;i<=n;i=x+1){x=(n/(n/i));ret=ad(ret,mod-1LL*(x-i+1)*Phi(n/i)%mod);}if(n<G){return p1[n]=ret;}else{return p2[m/n]=ret;}
}
map<int,int>mp[N];
int S(int n,int m){//cout<<" S "<if(m==0||n==0) return 0;if(n==1) return Phi(m);if(m==1) return phi[n];if(mp[n][m]) return mp[n][m];int ret=0;for(reg i=1;i*i<=n;++i){if(n%i==0){ret=ad(ret,(ll)phi[n/i]*S(i,m/i)%mod);if(i!=n/i) ret=ad(ret,(ll)phi[i]*S(n/i,m/(n/i))%mod);}}return mp[n][m]=ret;
}
int main(){rd(n);rd(m);if(n>m) swap(n,m);sieve(M-5);
//    cout<<" after sieve "<ll ans=0;memset(p1,-1,sizeof p1);memset(p2,-1,sizeof p2);for(reg i=1;i<=n;++i) {//    cout<<" ii "<ans=(ans+(ll)i/pro[i]*S(pro[i],m)%mod)%mod;}printf("%lld",ans);return 0;
}}
signed main(){Miracle::main();return 0;
}/*Author: *Miracle*Date: 2019/3/7 20:51:51
*/

转载于:https://www.cnblogs.com/Miracevin/p/10492891.html

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！