diff --git a/TangDou/Topic/SpareTire.cpp b/TangDou/Topic/SpareTire.cpp index 3e9fa04..4d0b6de 100644 --- a/TangDou/Topic/SpareTire.cpp +++ b/TangDou/Topic/SpareTire.cpp @@ -60,18 +60,19 @@ signed main() { t *= p[j]; } - // https://www.cnblogs.com/xiuwenli/p/9610029.html // 比如找到了s=6=2*3,需要知道s是奇数个,还是偶数个因子 // n/s:范围内6的倍数有多少个 int k = n / t; - int pt = k * (k + 1) % mod * (2 * k + 1) % mod * Six % mod; - pt = pt * t % mod * t % mod; - pt = (pt + k * (t + t * k) % mod * Two % mod) % mod; + int x = k * (k + 1) % mod * (2 * k + 1) % mod * Six % mod; + x = x * t % mod * t % mod; // 乘上t^2 + // 还需要累加等差数列部分 + // 首项是t,项数是k,末项 t*k + x = (x + k * (t + t * k) % mod * Two % mod) % mod; if (cnt & 1) - s = (s + pt) % mod; + s = (s + x) % mod; else - s = (s - pt + mod) % mod; + s = (s - x + mod) % mod; } // 输出 cout << (res - s + mod) % mod << endl; diff --git a/TangDou/Topic/容斥原理.md b/TangDou/Topic/容斥原理.md index e632bee..4c3256f 100644 --- a/TangDou/Topic/容斥原理.md +++ b/TangDou/Topic/容斥原理.md @@ -967,10 +967,11 @@ $a_{m+1}= (m+1)^2 + (m+1) $ 即$S_n=a_1+a_2+a_3+...+a_n \rightarrow$ $S_n=1^2+1+2^2+2+3^2+3+...+n^2+n$ 拆开来看: -① $1,2,3,...,n$很显然是一个等差数列,$S_n'=(1+n)*n/2$ +① $1,2,3,...,n$很显然是一个等差数列, +$$\large S_n'=(1+n)*n/2$$ ② $1^2,2^2,3^2,...,n^2$这个序列的求和公式是多少呢? 这是自然数平方求和数列,有公式: -$$\large S_n''=n*(n+1)(2*n+1)/6$$; +$$\large S_n''=n*(n+1)(2*n+1)/6$$ 这个东西是怎么来的呢? @@ -998,6 +999,21 @@ $S''=(n+1)(2n^2+n)/6$ 我们直接求与$m$互质的数较难,所以我们可以换个思路,求与 $m$不互质的数,那么与$m$不互质的数,是取$m$的素因子的乘积(因为根据唯一分解定理,任意个数都可看作的素数积),那么我们将$m$分解质因数,通过容斥定理,就可以得道与$m$不互质的数,总和$sum$减去这些数对应的$a$的和就是答案了。 +**代码细节** + +如果我们利用两个质数$2,3$组成了一个数$t=6$,那么在$1\sim n$范围内,一共几个$6$的倍数呢? 以前学习过,是 $n/6=n/t$ 个, +即:$t,2t,3t,…n/t*t$ + + +现在我们需要把$i \in [t,2t,3t,…n/t*t]$计算$\sum a[i]$ + +有公因子$t$,设$i=t*j$ +我们观察每一项: +$a(i)=i^2+i=(t*j)^2+(t*j)=t^2*j^2+t*j$ +在平方和公式前面,要乘一个系数$t$的平方,同时在等差数列求和公式前面要乘一个系数$t$。 + +根据通项公式,可以以$O(1)$时间快速计算出结果: +$$\large s_x=t^2*s_x'+t*s_x''$$ $Sample$ $Input$ ```cpp {.line-numbers} 4 4 @@ -1006,3 +1022,88 @@ $Sample$ $Output$ ```cpp {.line-numbers} 14 ``` + +#### $Code$ +```cpp {.line-numbers} +#include +using namespace std; +#define int long long +#define endl "\n" + +const int mod = 1e9 + 7; // 尽量这样定义mod ,减少非必要的麻烦 + +// 快速幂 +int qmi(int a, int b) { + int res = 1; + a %= mod; + while (b) { + if (b & 1) res = res * a % mod; + b >>= 1; + a = a * a % mod; + } + return res; +} + +vector p; // 将m拆分成的质数因子序列p + +signed main() { +#ifndef ONLINE_JUDGE + freopen("SpareTire.in", "r", stdin); +#endif + int n, m; + + int Six = qmi(6, mod - 2); // 因为需要用到 % 1e9+7 下6的逆元,用费马小定理+快速幂求逆元 + int Two = qmi(2, mod - 2); // 因为需要用到 % 1e9+7 下2的逆元,用费马小定理+快速幂求逆元 + + while (cin >> n >> m) { + // 所结果拆分为平方和公式,等差数列两部分 + // 注意:现在求的是整体值,还没有去掉不符合条件的数字 + int first = n * (n + 1) % mod * (2 * n + 1) % mod * Six % mod; + int second = n * (n + 1) % mod * Two % mod; + int res = (first + second) % mod; + + // 对m进行质因子分解 + int t = m; // 复制出来 + for (int i = 2; i * i <= t; i++) { + if (t % i == 0) { + p.push_back(i); + while (t % i == 0) t = t / i; + } + } + if (t > 1) p.push_back(t); + + /* + 容斥原理 + 例如有3个因子,那么1<<3=8(1000二进制) + 然后i从1开始枚举直到7(111二进制),i中二进制的位置1表式取这个位置的因子 + 例如i=3(11二进制) 表示取前两个因子,i=5(101)表示取第1个和第3个的因子 + */ + int s = 0; + for (int i = 1; i < (1 << p.size()); i++) { + int cnt = 0, t = 1; + for (int j = 0; j < p.size(); j++) + if ((i >> j) & 1) { + cnt++; + t *= p[j]; + } + + // 比如找到了s=6=2*3,需要知道s是奇数个,还是偶数个因子 + // n/s:范围内6的倍数有多少个 + int k = n / t; + int x = k * (k + 1) % mod * (2 * k + 1) % mod * Six % mod; + x = x * t % mod * t % mod; // 乘上t^2 + // 还需要累加等差数列部分 + // 首项是t,项数是k,末项 t*k + x = (x + k * (t + t * k) % mod * Two % mod) % mod; + + if (cnt & 1) + s = (s + x) % mod; + else + s = (s - x + mod) % mod; + } + // 输出 + cout << (res - s + mod) % mod << endl; + } +} + +``` \ No newline at end of file