'commit'

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;

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 <bits/stdc++.h>
+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<int> 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;
+    }
+}
+
+```