很显然我们可以推导出这个式子
设 \(a(m)\) 为 \(m\) 位置的值
\[ \mu(m) = \sum_{m | d} a(d) \a(m) = \sum_{m|d}\mu(\frac{d}{m})\mu(d) \a(m) = \sum_{i = 1}^{\lfloor \frac{n}{m} \rfloor} \mu(i)\mu(im) \a(m) = \mu(m) \sum_{i = 1}^{\lfloor \frac{n}{m} \rfloor} \mu(i)^{2}[gcd(m,i) == 1] \]
而 \(\mu(i)^{2}\) 的本质是无平方因子数, 这个可以容斥
容斥的方法是 (若没有其他限制)
\[ ans = \sum_{i = 1}^{\sqrt{N}} \mu(i)\lfloor \frac{N}{i^{2}}\rfloor \]
那么这里的就是
\[ ans = \sum_{i = 1}^{\sqrt{N / M}} [gcd(i,m) == 1]\mu(i)\sum_{j = 1}^{\lfloor \frac{N}{i ^ 2} \rfloor} [gcd(j,m) == 1] \]
前面的互质可以直接枚举
后面的互质可以通过莫比乌斯反演外加预处理 M 中莫比乌斯值不为 0 的数算出来
- #include <bits/stdc++.h>
- #define fi first
- #define se second
- #define pii pair<int,int>
- #define mp make_pair
- #define pb push_back
- #define space putchar(' ')
- #define enter putchar('\n')
- #define eps 1e-10
- #define ba 47
- #define MAXN 100005
- //#define ivorysi
- using namespace std;
- typedef long long int64;
- typedef unsigned int u32;
- typedef double db;
- template<class T>
- void read(T &res) {
- res = 0;T f = 1;char c = getchar();
- while(c <'0' || c> '9') {
- if(c == '-') f = -1;
- c = getchar();
- }
- while(c>= '0' && c <= '9') {
- res = res * 10 +c - '0';
- c = getchar();
- }
- res *= f;
- }
- template<class T>
- void out(T x) {
- if(x <0) {x = -x;putchar('-');}
- if(x>= 10) {
- out(x / 10);
- }
- putchar('0' + x % 10);
- }
- int64 N,MAXV;
- int M,mu[1000005];
- int prime[1000005],tot;
- bool nonprime[1000005];
- vector<pii> division;
- int gcd(int a,int b) {
- return b == 0 ? a : gcd(b,a % b);
- }
- int Mu(int x) {
- if(x <= 1000000) return mu[x];
- int res = 1;
- for(int i = 2 ; i <= x / i ; ++i) {
- if(x % i == 0) {
- int c = 0;
- while(x % i == 0) {x /= i;++c;}
- if(c>= 2) return 0;
- res = -res;
- }
- }
- if(x != 1) res = -res;
- return res;
- }
- int64 calc(int64 n) {
- int64 res = 0;
- for(auto t : division) {
- if(n <t.fi) break;
- res += 1LL * t.se * (n / t.fi);
- }
- return res;
- }
- void Init() {
- mu[1] = 1;
- for(int i = 2 ; i <= 1000000 ; ++i) {
- if(!nonprime[i]) {
- prime[++tot] = i;
- mu[i] = -1;
- }
- for(int j = 1 ; j <= tot ; ++j) {
- if(prime[j]> 1000000 / i) break;
- nonprime[i * prime[j]] = 1;
- if(i % prime[j] == 0) break;
- else mu[i * prime[j]] = -mu[i];
- }
- }
- }
- void Solve() {
- read(N);read(M);
- if(Mu(M) == 0) {puts("0");return;}
- division.clear();
- for(int i = 1 ; i <= M / i ; ++i) {
- if(M % i == 0) {
- int j = M / i;
- int x = Mu(i),y = Mu(j);
- if(x) division.pb(mp(i,x));
- if(i != j && y) division.pb(mp(j,y));
- }
- }
- sort(division.begin(),division.end());
- int64 T = N / M,res = 0;
- for(int i = 1 ; i <= T / i ; ++i) {
- if(gcd(i,M) == 1) {
- res += Mu(i) * calc(T / (i * i));
- }
- }
- res = res * Mu(M);
- out(res);enter;
- }
- int main(){
- #ifdef ivorysi
- freopen("f1.in","r",stdin);
- #endif
- Init();
- int T;
- read(T);
- for(int i = 1 ; i <= T ; ++i) Solve();
- }
来源: http://www.bubuko.com/infodetail-3090010.html