莫比乌斯反演部分习题

最近写反演题也快没头了，各种线性不会筛，各种卡常……

于是决定写一篇专题，来记录一下我写过的反演题目。

 

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\n",calc(a,b));
57     }
58     
59     return 0;
60 }
View Code

 

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\n",ans);
53     }
54     return 0;
55 }
View Code

 

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\n",calc(n,m));
59     }
60     
61     return 0;
62 }
View Code

 

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\n",ans[i]);
100     
101     return 0;
102 }
View Code

 

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\n",getans(n,m));
69     return 0;
70 }
View Code

 

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\n",query(n,m));
59     }
60     return 0;
61 }
View Code

 

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\n",query(a,b));
56     }
57     return 0;
58 }
View Code

 

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\n",sumphi(n),summu(n));
 98     }
 99     return 0;
100 }
View Code

 

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\n",calc(n,m));
65     }
66     return 0;
67 }
View Code

 

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\n",f(n));
76     return 0;
77 }
View Code

 

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\n",ans);
53      
54     return 0;
55 }
View Code

 

其实还有BZOJ3561，只不过弃坑大法好了。

关于杜教筛存储方式的正确性：你进去查询的时候总是令n/i=x的最小(或最大)的i，所以存储在位置n/i的方式正确。

posted @ 2018-01-07 21:10  Cmd2001  阅读(1714)  评论(0)  编辑  收藏  举报