TrickGCD（HDU 6053 莫比乌斯函数的反演）

TrickGCD

Time Limit: 5000/2500 MS (Java/Others)    Memory Limit: 262144/262144 K (Java/Others)

Total Submission(s): 2492    Accepted Submission(s): 958

Problem Description

You are given an array A , and Zhu wants to know there are how many different array B satisfy the following conditions?

* 1≤Bi≤Ai
* For each pair( l , r ) (1≤l≤r≤n) , gcd(bl,bl+1...br)≥2

 

Input

The first line is an integer T(1≤T≤10) describe the number of test cases.

Each test case begins with an integer number n describe the size of array A.

Then a line contains n numbers describe each element of A

You can assume that 1≤n,Ai≤105

 

Output

For the kth test case , first output "Case #k: " , then output an integer as answer in a single line . because the answer may be large , so you are only need to output answer mod 109+7

 

Sample Input

1
4
4 4 4 4

 

Sample Output

Case #1: 17

 

//题意：给你一个A数组，要求满足1<= B[i] <=A[i]且任意l，r，gcd(b[l],b[l+1],...,b[r])>=2的B数组的个数。

//思路：首先，我们要确定满足条件的gcd有k个，仔细想一下就知道是要求A数组里的最小值，2-min都可以，即

2<=k<=min。

然后，我们要算的是每个k的贡献值，val[k]=(A[1]/k)*(A[2]/k)*...*(A[n]/k)。

为什么是这样？其实很好理解，举个例子k=2，A[1]=4,A[2]=4,A[3]=4,A[4]=4，A[1]里可以取2或4都能满足k=2，A[i]里能有几个是可以满足k的？就是A[i]/k嘛。

但细心的小伙伴会发现这样是有问题的，上面那个例子val[2]里不是已经把val[4]的价值都包含进去了吗？是的，所有平方数的倍数val都等于0，因为都已经包含在前面了。

接下来思考要求的总数sum=所有贡献值之和吗？

显然不是，举个简单的例子，算k=2的贡献时把k=6时的贡献值已经算进去了，而算k=3的贡献的时候又把k=6的贡献算了一边，这样的重复情况显然还有很多，怎么办？这里用到的是容斥原理。

让我们先思考一下，什么样的数会被计算多次？

首先一个数可以分解成几个质数的乘积（就是因式分解嘛），如果一个数是有n个质数的成绩，那肯定要被计算多次了，每次算它的因数的时候都会被计算一次。

So，我们这样处理：

sum={ +val[k]，k是奇数个质数乘积；

            -val[k]，k是偶数个质数乘积；

           0，k是完全平方数的倍数；}

为什么奇数个质数乘积就要是+呢？举个例子，比如说30，30=2*3*5，在2的时候被 + 一次，在3的时候被 + 一次，在5的时候被 + 一次，6的时候被 - 一次，在10的时候被 - 一次，在15的时候被 - 一次，是不是相当于不加不减？！所以要 + 。

像4，9，12，16，18...这种完全平方数的倍数val就等于0。

这个+还是-还是0在coding的时候怎么实现？这里就用到了莫比乌斯函数。

Moblus(x)={  1，x=1；

                     (-1)^n，x=p1*p2*...*pn；（p为质数）

                     0，Rest；}

是不是觉得更上面提到的符号规律很像，就多了个负号（因为k=1的情况不考虑），这就是所谓的莫比乌斯函数的反演。（有模板）

以为这样就结束了吗？少年太naive！如果暴力算val的值是会TLE的！！！

这里要用筛法对k进行分块，用的筛法是埃氏筛法，即根据k的倍数来分块。

怎么实现？我们是设一个sum数组，一开始输入A[i]的时候，sum[A[i]]++，输入结束后再用个for循环（复杂度O（n）），是sum[i]里保存的是sum[0]—sum[i]的和。

此时，每个k的贡献val = 1^(sum[2*k-1]-sum[k-1])*2^(sum[3*k-1]-sum[2*k-1])*3^(sum[4*k-1]-sum[3*k-1])*...

为什么是这样呢？首先，如果gcd=k的话，那选的B[i]肯定是要比k大的（或等于），不然不可能有因子k。所以最小从sum[k-1]算起。而(k-1) -> (2*k-1) 的这段里是不能被k整除的（比k*1大，比k*2小），所以这段里都是只能*1的。同理，(2*k-1) -> (3*k-1) 这段是比k*2大，比k*3小的，都能选2个，所以都是*2的，即2的(sum[3*k-1]-sum[2*k-1])次方。之后也是同理。

计算a^b，在a，b都比较大的时候当然用快速幂啦！

这样就可以完美AC了！！！撒花！！！

PS：这是萌新看了网上代码一点点摸索出来的理解，有不对的地方欢迎指正。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <string>
#include <algorithm>
using namespace std;

const int MAX = 1e5 + 10;
const int mod = 1e9 + 7;

int n;
long long sum;
int a[MAX];
int ssum[MAX];
int Max, Min;

bool notp[MAX];
int prime[MAX], pnum, mu[MAX];

void Moblus()
{
	memset(notp, 0, sizeof(notp));
	mu[1] = 1;
	for (int i = 2; i < MAX; i++)
	{
		if (!notp[i])
		{
			prime[++pnum] = i;
			mu[i] = -1;
		}
		for (int j = 1; prime[j] * i < MAX; j++)
		{
			notp[prime[j] * i] = 1;
			if (i%prime[j] == 0)
			{
				mu[prime[j] * i] = 0;
				break;
			}
			mu[prime[j] * i] = -mu[i];
		}
	}
	for (int i = 0; i < MAX; i++)
		mu[i] = -mu[i];
}

long long fastpow(long long p, int m)
{
	if (p == 1 || m == 0)
		return 1;
	if (p == 0)
		return 0;
	long long res = 1;
	while (m)
	{
		if (m & 1)
			res = res*p%mod;
		p = p*p%mod;
		m >>= 1;
	}
	return res;
}

int main()
{
	int T, Case = 1;
	Moblus();
	scanf("%d", &T);
	while (T--)
	{
		sum = 0;
		scanf("%d", &n);
		memset(ssum, 0, sizeof(ssum));
		Min = MAX;
		Max = -1;
		for (int i = 0; i < n; i++)
		{
			scanf("%d", &a[i]);
			ssum[a[i]]++;
			if (Max < a[i])
				Max = a[i];
			if (a[i] < Min)
				Min = a[i];
		}
		printf("Case #%d: ", Case++);
		for (int i = 1; i <= Max; i++)
			ssum[i] += ssum[i - 1];

		long long res = 1;
		for (int i = 2; i <= Min; i++)
		{
			res = 1;
			if (mu[i] == 0)
				continue;
			int j = min(i, Max);
			int k = min(i * 2 - 1, Max);
			for (int p = 1;; p++)
			{
				if (ssum[k] - ssum[j - 1] != 0)
					res = res*fastpow(p, ssum[k] - ssum[j - 1]) % mod;
				if (k == Max)
					break;
				j += i;
				k += i;
				if (k > Max)
					k = Max;
			}
			sum = ((sum + mu[i] * res) % mod + mod) % mod;
		}

		printf("%lld\n", sum);
	}
	return 0;
}