python莫比乌斯内接矩形_莫比乌斯反演部分习题 - osc_9uwqh9yn的个人空间 - OSCHINA - 中文开源技术交流社区...

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

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

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

2 #include

3 #include

4 #include

5 #define lli long long int

6 using namespacestd;7 const int maxn=5e4+1e2;8

9 intmu[maxn];10 lli sum[maxn];11

12 inline voidgen() {13 static intprime[maxn],cnt;14 static unsigned charvis[maxn];15 sum[1] = mu[1] = 1;16 for(int i=2;i

33 inline lli calc(int n,intm) {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 returnret;43 }44

45 intmain() {46 static intT;47 inta,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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;8 const int maxn=5e4+1e2;9

10 lli summu[maxn],ans;11

12 inline void getmu(intlim) {13 static intprime[maxn],mu[maxn],cnt;14 static charvis[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 returnret;41 }42

43 intmain() {44 static intt,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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;8 const int maxn=1e7+1e2;9

10 lli f[maxn];11

12 inline voidgetf() {13 static intmu[maxn],prime[maxn],cnt;14 static charvis[maxn];15

16 mu[1] = 1;17 for(int i=2;i

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 returnret;48 }49

50 intmain() {51 static intt;52 staticlli 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

2 #include

3 #include

4 #include

5 #define debug cout

6 using namespacestd;7 const int maxn=1e5+1e2;8 const int mod=0x7fffffff;9

10 structBinaryIndexTree {11 intdat[maxn];12 #define lowbit(x) ( x & -x )

13 inline void update(int pos,intx) {14 while( pos

27 structQNode {28 intn,m,a,id;29 friend bool operator < (const QNode &a,const QNode &b) {30 return a.a

40 intmu[maxn];41 intans[maxn],n,m,t;42

43 inline voidpre() {44 static intprime[maxn],cnt;45 static charvis[maxn];46 mu[1] = 1;47 for(int i=2;i

67 inline int calc(int n,intm) {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 voidgetans() {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

90 intmain() {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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;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 intprime[maxn],cnt;16 static charvis[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 returnret;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 returnret;63 }64

65 intmain() {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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;8 const int maxn=1e7+1e2,mx=1e7;9 const int mod=1e8+9;10

11 lli sum[maxn];12

13 inline voidpre() {14 static intprime[maxn],cnt;15 static charvis[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 returnret;50 }51

52 intmain() {53 static intt,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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;8 const int maxn=1e7+1e2;9

10 lli sum[maxn];11

12 inline voidpre() {13 static intprime[maxn],most[maxn],pows[maxn],cnt;14 static charvis[maxn];15

16 for(lli i=2;i

19 prime[++cnt] =i;20 most[i] = 1 , pows[i] =i;21 }22 for(int j=1;j<=cnt&&i*prime[j]

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 returnret;47 }48

49 intmain() {50 pre();51 static intt,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

3 #define lli long long int

4 #define uint unsigned int

5 using namespacestd;6 const uint maxn=1677025,maxm=1295;7

8 lli phi[maxn],mu[maxn];9 lli memphi[maxm],memmu[maxm];10 unsigned charvisphi[maxn],vismu[maxm];11 uintt,n;12

13 inline void reset(unsigned char* dst,uintlen) {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 voidgen() {28 static uint prime[maxn/2],cnt;29 static unsigned charvis[maxn];30

31 phi[1] = mu[1] = 1;32 for(uint i=2;i

51 inline lli sumphi(uintx) {52 if( x > 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(uintx) {66 if( x

80 inline uintgetint() {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 returnret;88 }89

90 intmain() {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

2 #include

3 #include

4 #define lli long long int

5 #define debug cout

6 using namespacestd;7 const int maxn=5e6+1e2;8 const int mod=1e9+7;9

10 lli sum[maxn];11 intk;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 returnret;23 }24 inline voidinit() {25 static intpows[maxn],prime[maxn],cnt;26 static charvis[maxn];27 pows[1] = sum[1] = 1;28 for(lli i=2;i

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 returnret;56 }57

58 intmain() {59 static intT,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

2 #include

3 #include

4 #include

5 #define lli long long int

6 #define debug cout

7 using namespacestd;8 const int maxn=1e6+1e2,lim=1e6;9 const int mod = 1000000007;10

11 lli mu[maxn],mem[maxn];12 intn;13

14 inline voidpre() {15 static intprime[maxn],cnt;16 static charvis[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(intp) {37 return p<=lim ? mu[p] : mem[n/p];38 }39 inline voidsieve() {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 returnret;69 }70

71 intmain() {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

2 #include

3 #include

4 #include

5 #define lli long long int

6 using namespacestd;7 const int maxn=1e5+1e2;8

9 lli phi[maxn];10 lli ans;11 intn,m;12

13 inline voidgen()14 {15 static intprime[maxn],cnt;16 static unsigned charvis[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 intmain()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的方式正确。