这题很容易转化到一个容斥计数问题。而用指数复杂度的枚举计数法显然会挂,只能考虑别的方法。
首先将a[i]用gcd(a[i], m)替换,排序去重后得到一组m的约数,而m不超过1e9,因此m的所有约数最多在1000左右。
假设数组a只含有2,3两个元素,那么显然答案应该是f(2) + f(3) - f(6),f(i)表示只有i的情况下的答案。
2和3的系数是1,6的系数是-1,其余系数为0。
根据容斥原理:对于某个特定的约数,若存在数组a的一个子集其lcm为该约数,那么若集合元素个数为偶数,贡献-1,若为奇数,贡献1。
而最终的答案一定由所有约数合成,其系数取值{-1,0,1}。
因此我们考虑数组a的前m个数,其对应的系数分布是可以确定的。
加入第m+1个数时,显然新的集合必然由原集合的子集并上新元素得到,分别考虑m的所有约数k,其对应的系数为op(k)
若op(k)为-1,则说明偶数子集的个数比奇数子集的个数多1,若我们添如新元素x,则恰好使得奇数子集个数比偶数多1,但此时贡献
不是对k的,而是对lcm(k, x),用公式表示就是:
op(lcm(k, x)) += op(k) * (-1).
其余情形类似。
最终答案是ans = sigma(op(i) * f(i)), i | m.
这样我们只需最多进行1000次遍历就可使得元素个数加一,总复杂度约为O(1e6)。
1 #include <cstdio> 2 #include <cstring> 3 #include <algorithm> 4 #include <map> 5 #include <string> 6 #include <vector> 7 #include <set> 8 #include <cmath> 9 #include <ctime> 10 #include <cassert> 11 #define lson (u << 1) 12 #define rson (u << 1 | 1) 13 #define cls(i, j) memset(i, j, sizeof i) 14 #pragma comment(linker, "/STACK:102400000,102400000") 15 using namespace std; 16 typedef __int64 ll; 17 const double eps = 1e-6; 18 const double pi = acos(-1.0); 19 const int maxn = 1e4 + 10; 20 const int maxm = 4e4 + 10; 21 const int inf = 0x3f3f3f3f; 22 const ll linf = 0x3fffffffffffffff; 23 const ll mod = 1e9 + 7; 24 25 map<ll, int> mapi; 26 int n; 27 ll a[maxn], m; 28 int prime[maxm], k; 29 bool vis[maxm]; 30 ll fact[maxn], nf; 31 int cnt[maxn]; 32 ll table[maxn], nt; 33 int op[maxn]; 34 int op1[maxn]; 35 ll ans; 36 37 ll gcd(ll a, ll b) { return !b ? a : gcd(b, a % b); } 38 39 void init(){ 40 cls(vis, 0); 41 k = 0; 42 for(int i = 2; i < maxm; i++){ 43 if(vis[i]) continue; 44 prime[k++] = i; 45 for(int j = i * 2; j < maxm; j += i) vis[j] = 1; 46 } 47 } 48 49 void dfs(int next, ll num){ 50 if(next >= nf){ 51 table[nt++] = num; 52 mapi[num] = nt - 1; 53 return; 54 } 55 ll tem = num; 56 for(int i = 0; i <= cnt[next]; i++){ 57 dfs(next + 1, tem); 58 tem *= fact[next]; 59 } 60 } 61 62 ll f(ll num){ 63 ll tem = (m - 1) / num; 64 return tem * (tem + 1) / 2 * num; 65 } 66 67 void factorize(){ 68 nf = 0; 69 ll m1 = m; 70 int mid = (int)sqrt(m); 71 for(int i = 0; prime[i] <= mid; i++){ 72 if(m % prime[i]) continue; 73 fact[nf++] = prime[i]; 74 cnt[nf - 1] = 0; 75 while(m % prime[i] == 0) ++cnt[nf - 1], m /= prime[i]; 76 mid = (int)sqrt(m); 77 } 78 if(m != 1) fact[nf++] = m, cnt[nf - 1] = 1; 79 m = m1; 80 } 81 82 ll solve(){ 83 mapi.clear(); 84 int n1 = 0; 85 for(int i = 0; i < n; i++){ 86 ll tem = a[i] % m; 87 if(!tem) continue; 88 a[n1++] = gcd(m, tem); 89 } 90 n = n1; 91 sort(a, a + n); 92 n = unique(a, a + n) - a; 93 cls(vis, 0); 94 for(int i = 0; i < n; i++){ 95 for(int j = i + 1; j < n; j++){ 96 if(a[i] % a[j] == 0) vis[i] = 1; 97 else if(a[j] % a[i] == 0) vis[j] = 1; 98 } 99 } 100 n1 = 0; 101 for(int i = 0; i < n; i++) if(!vis[i]) a[n1++] = a[i]; 102 n = n1; 103 sort(a, a + n); 104 factorize(); 105 nt = 0; 106 dfs(0, 1); 107 cls(op, 0); 108 ans = 0; 109 for(int i = 0; i < n; i++){ 110 cls(op1, 0); 111 for(int j = 0; j < nt; j++){ 112 ll tem = a[i] / gcd(a[i], table[j]) * table[j]; 113 if(tem >= m) continue; 114 op1[mapi[tem]] += -op[j]; 115 } 116 ++op1[mapi[a[i]]]; 117 for(int j = 0; j < nt; j++) op[j] += op1[j]; 118 } 119 // for(int i = 0; i < nt; i++) printf("%d ", op[i]); 120 // puts(""); 121 for(int i = 0; i < nt; i++) ans += op[i] * f(table[i]); 122 return ans; 123 } 124 125 int main(){ 126 // freopen("in.txt", "r", stdin); 127 int T, kase = 0; 128 init(); 129 scanf("%d", &T); 130 while(T--){ 131 scanf("%d%I64d", &n, &m); 132 for(int i = 0; i < n; i++) scanf("%I64d", &a[i]); 133 ll ans = solve(); 134 printf("Case #%d: %I64d ", ++kase, ans); 135 } 136 return 0; 137 }