最近写反演题也快没头了，各种线性不会筛，各种卡常……
于是决定写一篇专题，来记录一下我写过的反演题目。
BZOJ1101：
求1<=i<=n,1<=j<=m，gcd(i,j)==d的对数。
先让n/=d,m/=d，变成求gcd(i,j)==1的对数。
然后预处理出μ(d)的前缀和，O(sqrt(n))枚举d即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 using namespace std; 7 const int maxn=5e4+1e2; 8 9 int mu[maxn]; 10 lli sum[maxn]; 11 12 inline void gen() { 13 static int prime[maxn],cnt; 14 static unsigned char vis[maxn]; 15 sum[1] = mu[1] = 1; 16 for(int i=2;i<maxn;i++) { 17 if( !vis[i] ) { 18 prime[++cnt] = i, 19 mu[i] = -1; 20 } 21 for(int j=1;j<=cnt&&i*prime[j]<maxn;j++) { 22 vis[i*prime[j]] = 1; 23 mu[i*prime[j]] = -mu[i]; 24 if( ! ( i % prime[j]) ) { 25 mu[i*prime[j]] = 0; 26 break; 27 } 28 } 29 sum[i] = sum[i-1] + mu[i]; 30 } 31 } 32 33 inline lli calc(int n,int m) { 34 lli ret = 0; 35 if( n > m ) 36 swap(n,m); 37 int pos = 0; 38 for(int i=1;i<=n;i=pos+1) { 39 pos = min( n / ( n / i ) , m / ( m / i ) ); 40 ret += ( sum[pos] - sum[i-1] ) * ( n / i ) * ( m / i ); 41 } 42 return ret; 43 } 44 45 int main() { 46 static int T; 47 int a,b,d; 48 49 gen(); 50 51 scanf("%d",&T); 52 53 while( T-- ) { 54 scanf("%d%d%d",&a,&b,&d); 55 a /= d , b /= d; 56 printf("%lld ",calc(a,b)); 57 } 58 59 return 0; 60 }
BZOJ2301：
求二维区间gcd(a,b)==k，利用二维前缀和思想转化为1101即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=5e4+1e2; 9 10 lli summu[maxn],ans; 11 12 inline void getmu(int lim) { 13 static int prime[maxn],mu[maxn],cnt; 14 static char vis[maxn]; 15 mu[1] = 1; 16 for(int i=2;i<=lim;i++) { 17 if( !vis[i] ) 18 prime[++cnt] = i, 19 mu[i] = -1; 20 for(int j=1;j<=cnt&&i*prime[j]<=lim;j++) { 21 vis[i*prime[j]] = 1; 22 mu[i*prime[j]] = -mu[i]; 23 if( ! ( i % prime[j] ) ) { 24 mu[i*prime[j]] = 0; 25 break; 26 } 27 } 28 } 29 for(int i=1;i<=lim;i++) 30 summu[i] = summu[i-1] + mu[i]; 31 } 32 33 inline lli calc(lli n,lli m) { 34 lli ret = 0; 35 const lli t = min(n,m); 36 for(lli i=1,j;i<=t;i=j+1) { 37 j = min( n / ( n / i ) , m / ( m / i ) ); 38 ret += ( summu[j] - summu[i-1] ) * ( n / i ) * ( m / i ); 39 } 40 return ret; 41 } 42 43 int main() { 44 static int t,a,b,c,d,k; 45 getmu(5e4); 46 scanf("%d",&t); 47 while( t-- ) { 48 scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); 49 a-- , c--; 50 a /= k , b /= k , c /= k , d /= k; 51 ans = calc(b,d) + calc(a,c) - calc(b,c) - calc(a,d); 52 printf("%lld ",ans); 53 } 54 return 0; 55 }
BZOJ2820：
求二维gcd为质数。是质数这一点不太好处理，暂且留着它。
然后我们nlogn筛出后面那个东西，跑前缀和，O(sqrt(n))枚举p即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=1e7+1e2; 9 10 lli f[maxn]; 11 12 inline void getf() { 13 static int mu[maxn],prime[maxn],cnt; 14 static char vis[maxn]; 15 16 mu[1] = 1; 17 for(int i=2;i<maxn;i++) { 18 if( !vis[i] ) { 19 prime[++cnt] = i, 20 mu[i] = -1; 21 } 22 for(int j=1;j<=cnt&&i*prime[j]<maxn;j++) { 23 vis[i*prime[j]] = 1; 24 mu[i*prime[j]] = -mu[i]; 25 if( !( i % prime[j] ) ) { 26 mu[i*prime[j]] = 0; 27 break; 28 } 29 } 30 } 31 for(int i=1;i<=cnt;i++) { 32 const int p = prime[i]; 33 for(int j=1;j*p<maxn;j++) 34 f[j*p] += mu[j]; 35 } 36 for(int i=1;i<maxn;i++) 37 f[i] += f[i-1]; 38 } 39 40 inline lli calc(lli n,lli m) { 41 const lli mi = min(n,m); 42 lli ret = 0; 43 for(lli i=1,j;i<=mi;i=j+1) { 44 j = min( n / ( n / i ) , m / ( m / i ) ); 45 ret += ( f[j] - f[i-1] ) * ( n / i ) * ( m / i ); 46 } 47 return ret; 48 } 49 50 int main() { 51 static int t; 52 static lli n,m; 53 getf(); 54 55 scanf("%d",&t); 56 while( t-- ) { 57 scanf("%lld%lld",&n,&m); 58 printf("%lld ",calc(n,m)); 59 } 60 61 return 0; 62 }
BZOJ3529：
我们能枚举gcd，算出gcd为(i)的组数g(i)为：
对于a没有限制，我们能化简为：
注意到n,m范围很小，我们能nlogn筛出后面的那个卷积，然后按照a递增插入坐标为d的树状数组并用树状数组的前缀和计算答案。复杂度nlog^2n。因为此题卡常，所以用int自然溢出取模。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define debug cout 6 using namespace std; 7 const int maxn=1e5+1e2; 8 const int mod=0x7fffffff; 9 10 struct BinaryIndexTree { 11 int dat[maxn]; 12 #define lowbit(x) ( x & -x ) 13 inline void update(int pos,int x) { 14 while( pos < maxn ) 15 dat[pos] += x, 16 pos += lowbit(pos); 17 } 18 inline int query(int pos) { 19 int ret = 0; 20 while( pos ) 21 ret += dat[pos], 22 pos -= lowbit(pos); 23 return ret; 24 } 25 }bit; 26 27 struct QNode { 28 int n,m,a,id; 29 friend bool operator < (const QNode &a,const QNode &b) { 30 return a.a < b.a; 31 } 32 }ns[maxn]; 33 struct FNode { 34 int f,div; 35 friend bool operator < (const FNode &a,const FNode &b) { 36 return a.f < b.f; 37 } 38 }fs[maxn]; 39 40 int mu[maxn]; 41 int ans[maxn],n,m,t; 42 43 inline void pre() { 44 static int prime[maxn],cnt; 45 static char vis[maxn]; 46 mu[1] = 1; 47 for(int i=2;i<maxn;i++) { 48 if( !vis[i] ) 49 prime[++cnt] = i, 50 mu[i] = -1; 51 for(int j=1;j<=cnt&&i*prime[j]<maxn;j++) { 52 vis[i*prime[j]] = 1, 53 mu[i*prime[j]] = -mu[i]; 54 if( ! ( i % prime[j] ) ) { 55 mu[i*prime[j]] = 0; 56 break; 57 } 58 } 59 } 60 for(int i=1;i<maxn;i++) { 61 fs[i].div = i; 62 for(int j=i;j<maxn;j+=i) 63 fs[j].f += i; 64 } 65 } 66 67 inline int calc(int n,int m) { 68 int ret = 0 , lim = min(n,m); 69 for(int i=1,j;i<=lim;i=j+1) { 70 j = min( n / ( n / i ) , m / ( m / i ) ); 71 ret += ( n / i ) * ( m / i ) * ( bit.query(j) - bit.query(i-1) ); 72 } 73 return ret & mod; 74 } 75 76 inline void getans() { 77 sort(ns+1,ns+1+n); 78 sort(fs+1,fs+maxn); 79 int pos = 0; 80 for(int i=1;i<=n;i++) { 81 while( pos + 1 < maxn && fs[pos+1].f <= ns[i].a ) { 82 ++pos; 83 for(int j=fs[pos].div;j<maxn;j+=fs[pos].div) 84 bit.update(j,fs[pos].f*mu[j/fs[pos].div]); 85 } 86 ans[ns[i].id] = calc(ns[i].n,ns[i].m); 87 } 88 } 89 90 int main() { 91 scanf("%d",&n); 92 for(int i=1;i<=n;i++) 93 scanf("%d%d%d",&ns[i].n,&ns[i].m,&ns[i].a), 94 ns[i].id = i; 95 pre(); 96 getans(); 97 98 for(int i=1;i<=n;i++) 99 printf("%d ",ans[i]); 100 101 return 0; 102 }
BZOJ2154：
然后后面等差数列求和，筛出μ(d)*d*d即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=1e7+1e2; 9 const int mod = 20101009; 10 11 lli sum[maxn],ss[maxn]; 12 lli n,m; 13 14 inline void pre(lli lim) { // don't prepare to maxn. 15 static int prime[maxn],cnt; 16 static char vis[maxn]; 17 18 sum[1] = 1; 19 for(lli i=2;i<=lim;i++) { 20 if( !vis[i] ) { 21 prime[++cnt] = i, 22 sum[i] = -1; 23 } 24 for(int j=1;j<=cnt&&i*prime[j]<=lim;j++) { 25 vis[i*prime[j]] = 1; 26 sum[i*prime[j]] = -sum[i]; 27 if( ! ( i % prime[j]) ) { 28 sum[i*prime[j]] = 0; 29 break; 30 } 31 } 32 } 33 for(int i=1;i<=lim;i++) 34 sum[i] = sum[i] * i * i % mod; 35 for(int i=1;i<=lim;i++) { 36 sum[i] += sum[i-1], 37 sum[i] = ( sum[i] % mod + mod ) % mod; 38 ss[i] = ( ss[i-1] + i ) % mod; 39 } 40 } 41 42 inline lli fnm(lli n,lli m,lli step) { 43 n /= step , m /= step; 44 return ( ( n + 1 ) * n >> 1 ) % mod * ( ( ( m + 1 ) * m >> 1 ) % mod ) % mod; 45 } 46 inline lli f(lli n,lli m) { 47 lli ret = 0 , lim = min( n , m ); 48 for(lli i=1,j;i<=lim;i=j+1) { 49 j = min( n / ( n / i ) , m / ( m / i ) ); 50 ret += ( ( sum[j] - sum[i-1] ) % mod + mod ) % mod * fnm(n,m,i) % mod; 51 ret = ( ret % mod + mod ) % mod; 52 } 53 return ret; 54 } 55 inline lli getans(lli n,lli m) { 56 lli ret = 0 , lim = min( n , m ); 57 for(lli i=1,j;i<=lim;i=j+1) { 58 j = min( n / ( n / i ) , m / ( m / i ) ); 59 ret += ( ( ss[j] - ss[i-1] ) % mod + mod ) % mod * f(n/i,m/i); 60 ret = ( ret % mod + mod ) % mod; 61 } 62 return ret; 63 } 64 65 int main() { 66 scanf("%lld%lld",&n,&m); 67 pre(min(n,m)); 68 printf("%lld ",getans(n,m)); 69 return 0; 70 }
BZOJ2693：
同上一题公式，多次询问。
我们定义
然后这东西是这样子的：
然后后面的东西是可以线性筛的，当i%prime[j]不为0时直接相乘，等于0是f[i*prime[j]]=f[i]，因为新加入的因子全被μ消掉了，剩下的还是原来的因子，自己画一下就明白了QAQ。
其实不用这样分析的，写一个暴力筛看一看，基本上前10个数就能把WA的筛法全都WA掉。实在不行交nlogn的也能拿到不少分数。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=1e7+1e2,mx=1e7; 9 const int mod=1e8+9; 10 11 lli sum[maxn]; 12 13 inline void pre() { 14 static int prime[maxn],cnt; 15 static char vis[maxn]; 16 sum[1] = 1; 17 for(lli i=2;i<=mx;i++) { 18 if( !vis[i] ) { 19 prime[++cnt] = i; 20 sum[i] = 1 - i; 21 } 22 for(int j=1;j<=cnt&&i*prime[j]<=mx;j++) { 23 vis[i*prime[j]] = 1; 24 sum[i*prime[j]] = sum[i] * sum[prime[j]] % mod; 25 if( ! ( i % prime[j]) ) { 26 sum[i*prime[j]] = ( sum[i] ) % mod; 27 break; 28 } 29 } 30 } 31 for(int i=1;i<=mx;i++) { 32 sum[i] = sum[i] * i % mod, 33 sum[i] += sum[i-1], 34 sum[i] = ( sum[i] % mod + mod ) % mod; 35 } 36 } 37 38 inline lli calc(lli n,lli m,lli step) { 39 n /= step , m /= step; 40 return ( ( ( n + 1 ) * n >> 1 ) % mod ) * ( ( ( m + 1 ) * m >> 1 ) % mod ) % mod; 41 } 42 inline lli query(lli n,lli m) { 43 lli ret = 0 , lim = min( n , m ); 44 for(lli i=1,j;i<=lim;i=j+1) { 45 j = min( n / ( n / i ) , m / ( m / i ) ); 46 ret += calc(n,m,i) * ( ( sum[j] - sum[i-1] ) % mod + mod ) % mod; 47 ret = ( ret % mod + mod ) % mod; 48 } 49 return ret; 50 } 51 52 int main() { 53 static int t,n,m; 54 pre(); 55 scanf("%d",&t); 56 while( t-- ) { 57 scanf("%d%d",&n,&m); 58 printf("%lld ",query(n,m)); 59 } 60 return 0; 61 }
BZOJ3309：
右边的东西F(g)当然还是线性筛，分析一发发现：
对于g的唯一分解，如果所有质数的次数不相同，那么右面的F(g)为0，因为我们可以通过改变最小质数的位置使之为0，如果不全部相同，那么就会有一个是-a，一个是a-1(这时d=所有质数次数-1的连乘)的情况，这样消出来是-1。考虑上μ的影响，数值为(-1)^(k+1)，k为分解后的质数个数。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=1e7+1e2; 9 10 lli sum[maxn]; 11 12 inline void pre() { 13 static int prime[maxn],most[maxn],pows[maxn],cnt; 14 static char vis[maxn]; 15 16 for(lli i=2;i<maxn;i++) { 17 if( !vis[i] ) { 18 sum[i] = 1; // if i is a prime then miu(1) = 1 and f(i) = 1 , miu(i) = -1 but f(1) = 0 , so it should be 1. 19 prime[++cnt] = i; 20 most[i] = 1 , pows[i] = i; 21 } 22 for(int j=1;j<=cnt&&i*prime[j]<maxn;j++) { 23 const int tar = i * prime[j]; 24 vis[tar] = 1; 25 most[tar] = 1 , pows[tar] = prime[j]; 26 sum[tar] = most[i] == 1 ? -sum[i] : 0; 27 if( ! ( i % prime[j]) ) { 28 most[tar] = most[i] + 1, 29 pows[tar] = pows[i] * prime[j]; 30 sum[tar] = tar==pows[tar] ? 1 : most[tar/pows[tar]] == most[tar] ? -sum[tar/pows[tar]] : 0; 31 break; 32 } 33 } 34 } 35 for(int i=1;i<maxn;i++) 36 sum[i] += sum[i-1]; 37 } 38 39 inline lli query(lli n,lli m) { 40 const lli lim = min(n,m); 41 lli ret = 0; 42 for(lli i=1,j;i<=lim;i=j+1) { 43 j = min( n / ( n / i ) , m / ( m / i ) ); 44 ret += ( sum[j] - sum[i-1] ) * ( n / i ) * ( m / i ); 45 } 46 return ret; 47 } 48 49 int main() { 50 pre(); 51 static int t,a,b; 52 scanf("%d",&t); 53 while(t--) { 54 scanf("%d%d",&a,&b); 55 printf("%lld ",query(a,b)); 56 } 57 return 0; 58 }
BZOJ3944：
杜教筛裸题，自行学习吧。
另外，此题卡常，不要用map存。我是用数组+vis标记存的。清空的时候可以用达夫机器卡一下常，比memset快2333。
代码：
1 #pragma GCC optimize(2) 2 #include<cstdio> 3 #define lli long long int 4 #define uint unsigned int 5 using namespace std; 6 const uint maxn=1677025,maxm=1295; 7 8 lli phi[maxn],mu[maxn]; 9 lli memphi[maxm],memmu[maxm]; 10 unsigned char visphi[maxn],vismu[maxm]; 11 uint t,n; 12 13 inline void reset(unsigned char* dst,uint len) { 14 uint ite = ( len + 7 ) >> 3; 15 switch( len & 7u ) { 16 case 0 : do{ *dst++ = 0; 17 case 7 : *dst++ = 0; 18 case 6 : *dst++ = 0; 19 case 5 : *dst++ = 0; 20 case 4 : *dst++ = 0; 21 case 3 : *dst++ = 0; 22 case 2 : *dst++ = 0; 23 case 1 : *dst++ = 0; }while(--ite); 24 } 25 } 26 27 inline void gen() { 28 static uint prime[maxn/2],cnt; 29 static unsigned char vis[maxn]; 30 31 phi[1] = mu[1] = 1; 32 for(uint i=2;i<maxn;i++) { 33 if( !vis[i] ) { 34 prime[++cnt] = i, 35 phi[i] = i - 1 , mu[i] = -1; 36 } 37 for(uint j=1;j<=cnt&&((unsigned lli)i)*prime[j]<maxn;j++) { 38 const uint tar = i*prime[j]; 39 vis[tar] = 1; 40 if( ! ( i % prime[j] ) ) { 41 phi[tar] = phi[i] * prime[j] , mu[tar] = 0; 42 break; 43 } 44 phi[tar] = phi[i] * ( prime[j] - 1 ) , mu[tar] = -mu[i]; 45 } 46 } 47 for(uint i=1;i<maxn;i++) 48 phi[i] += phi[i-1] , mu[i] += mu[i-1]; 49 } 50 51 inline lli sumphi(uint x) { 52 if( x < maxn ) 53 return phi[x]; 54 uint mp = n / x; 55 if( visphi[mp] ) 56 return memphi[mp]; 57 lli ret = ( (lli) x ) * ( x + 1 ) >> 1; 58 for(uint i=2,j;i<=x;i=j+1) { 59 j = x / ( x / i ); 60 ret -= ( (lli) ( j - i + 1 ) ) * sumphi( x / i ); 61 } 62 visphi[mp] = 1; 63 return memphi[mp] = ret; 64 } 65 inline lli summu(uint x) { 66 if( x < maxn ) 67 return mu[x]; 68 uint mp = n / x; 69 if( vismu[mp] ) 70 return memmu[mp]; 71 lli ret = 1; 72 for(uint i=2,j;i<=x;i=j+1) { 73 j = x / ( x / i ); 74 ret -= ( (lli) ( j - i + 1 ) ) * summu( x / i ); 75 } 76 vismu[mp] = 1; 77 return memmu[mp] = ret; 78 } 79 80 inline uint getint() { 81 uint ret = 0 , ch = getchar(); 82 while( ch < '0' || ch > '9' ) 83 ch = getchar(); 84 while( '0' <= ch && ch <= '9' ) 85 ret = ret * 10 + ( ch - '0' ), 86 ch = getchar(); 87 return ret; 88 } 89 90 int main() { 91 gen(); 92 t = getint(); 93 while( t-- ) { 94 reset(visphi,maxm); 95 reset(vismu,maxm); 96 n = getint(); 97 printf("%lld %lld ",sumphi(n),summu(n)); 98 } 99 return 0; 100 }
BZOJ4407：
注意题目的k是不变的，否则不可做了……
日常反演，成：
(什么你推不出来这一步？回去做前面的题吧孩子)
然后考虑右边怎么筛，可以消去一个d然后用我在3309用的筛法，或者直接筛。
当i%prime[j]不为0时直接相乘，等于0是f[i*prime[j]]=f[i]*prime[j]^k。
为什么？考虑μ，新的含p^2的因数全为0了。其余的一些可以消掉的说。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<algorithm> 4 #define lli long long int 5 #define debug cout 6 using namespace std; 7 const int maxn=5e6+1e2; 8 const int mod=1e9+7; 9 10 lli sum[maxn]; 11 int k; 12 13 inline lli fastpow(lli base) { 14 int tme = k; 15 lli ret = 1; 16 while( tme ) { 17 if( tme & 1 ) 18 ret = ret * base % mod; 19 base = base * base % mod; 20 tme >>= 1; 21 } 22 return ret; 23 } 24 inline void init() { 25 static int pows[maxn],prime[maxn],cnt; 26 static char vis[maxn]; 27 pows[1] = sum[1] = 1; 28 for(lli i=2;i<maxn;i++) { 29 if( !vis[i] ) { 30 prime[++cnt] = i; 31 pows[i] = fastpow(i); 32 sum[i] = pows[i] - 1; 33 } 34 for(int j=1;j<=maxn&&i*prime[j]<maxn;j++) { 35 const int tar = i * prime[j]; 36 vis[tar] = 1 ; 37 if( ! ( i % prime[j] ) ) { 38 sum[tar] = sum[i] * pows[prime[j]] % mod; 39 break; 40 } 41 sum[tar] = sum[i] * sum[prime[j]] % mod; 42 } 43 } 44 for(int i=1;i<maxn;i++) 45 sum[i] += sum[i-1]; 46 } 47 48 inline lli calc(lli n,lli m) { 49 lli ret = 0 , lim = min(n,m); 50 for(lli i=1,j;i<=lim;i=j+1) { 51 j = min( n / ( n / i ) , m / ( m / i ) ); 52 ret += ( ( sum[j] - sum[i-1] ) % mod + mod ) % mod * ( n / i ) % mod * ( m / i ) % mod; 53 ret = ( ret % mod + mod ) % mod; 54 } 55 return ret; 56 } 57 58 int main() { 59 static int T,n,m; 60 scanf("%d%d",&T,&k); 61 init(); 62 while( T-- ) { 63 scanf("%d%d",&n,&m); 64 printf("%lld ",calc(n,m)); 65 } 66 return 0; 67 }
BZOJ4176：
杜教筛μ即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 #define debug cout 7 using namespace std; 8 const int maxn=1e6+1e2,lim=1e6; 9 const int mod = 1000000007; 10 11 lli mu[maxn],mem[maxn]; 12 int n; 13 14 inline void pre() { 15 static int prime[maxn],cnt; 16 static char vis[maxn]; 17 mu[1] = 1; 18 for(lli i=2;i<=lim;i++) { 19 if( !vis[i] ) { 20 prime[++cnt] = i; 21 mu[i] = -1; 22 } 23 for(int j=1;j<=cnt&&i*prime[j]<=lim;j++) { 24 const int tar = i * prime[j]; 25 vis[tar] = 1; 26 if( ! ( i % prime[j] ) ) { 27 mu[tar] = 0; 28 break; 29 } 30 mu[tar] = -mu[i]; 31 } 32 } 33 for(int i=1;i<=lim;i++) 34 mu[i] += mu[i-1]; 35 } 36 inline lli getmu(int p) { 37 return p<=lim ? mu[p] : mem[n/p]; 38 } 39 inline void sieve() { 40 int t = ( n + lim - 1 ) / lim; 41 while( t ) { 42 const int m = n / t; 43 mem[t] = 1; 44 for(lli i=2,j;i<=m;i=j+1) { 45 j = m / ( m / i ); 46 mem[t] -= ( j - i + 1 ) * ( getmu(m/j) ) % mod; 47 } 48 mem[t] %= mod; 49 t--; 50 } 51 } 52 inline lli sum(lli n) { 53 lli ret = 0; 54 for(lli i=1,j;i<=n;i=j+1) { 55 j = n / ( n / i ); 56 ret += ( j - i + 1 ) * ( n / i ) % mod, 57 ret %= mod; 58 } 59 return ret * ret % mod; 60 } 61 inline lli f(lli n) { 62 lli ret = 0; 63 for(lli i=1,j;i<=n;i=j+1) { 64 j = n / ( n / i ); 65 ret += ( getmu(j) - getmu(i-1) ) * sum( n / i ) % mod; 66 ret = ( ret % mod + mod ) % mod; 67 } 68 return ret; 69 } 70 71 int main() { 72 scanf("%d",&n); 73 pre(); 74 sieve(); 75 printf("%lld ",f(n)); 76 return 0; 77 }
BZOJ2005：
公式为：
只考虑含gcd的部分就可以了。
把gcd变成id，根据id=phi*1进行反演。
然后线性筛φ即可。
代码：
1 #include<iostream> 2 #include<cstdio> 3 #include<cstring> 4 #include<algorithm> 5 #define lli long long int 6 using namespace std; 7 const int maxn=1e5+1e2; 8 9 lli phi[maxn]; 10 lli ans; 11 int n,m; 12 13 inline void gen() 14 { 15 static int prime[maxn],cnt; 16 static unsigned char vis[maxn]; 17 18 if( n > m ) // m is the limit number 19 swap(n,m); 20 21 phi[1] = 1; 22 23 for(int i=2;i<=m;i++) 24 { 25 if( !vis[i] ) 26 prime[++cnt] = i, 27 phi[i] = i-1; 28 for(int j=1;j<=cnt&&i*prime[j]<=m;j++) 29 { 30 vis[i*prime[j]] = 1; 31 phi[i*prime[j]] = phi[i] * ( prime[j] - 1 ); 32 if( ! ( i % prime[j]) ) 33 { 34 phi[i*prime[j]] = phi[i] * prime[j]; 35 break; 36 } 37 } 38 } 39 } 40 41 int main() 42 { 43 scanf("%d%d",&n,&m); 44 45 gen(); 46 47 for(int i=1;i<=n;i++) 48 ans += phi[i] * ( n / i ) * ( m / i ); 49 50 ans = ans * 2 - ( (long long) n ) * m; 51 52 printf("%lld ",ans); 53 54 return 0; 55 }
其实还有BZOJ3561，只不过弃坑大法好了。
关于杜教筛存储方式的正确性：你进去查询的时候总是令n/i=x的最小(或最大)的i，所以存储在位置n/i的方式正确。