利用容斥原理求解 [a,b] 区间中与 n 互质的数字个数

文章目录

• 题目描述
• 分析
• 代码
• 拓展：ZOJ 3547

先看这道丧心病狂的题目：HDU 4135 Co-prime。题目大意就是，有 $T$ 组询问，每组询问给你三个数：$a, b, c$，问你闭区间 $[a,b]$ 中有多少个数字与 $n$ 互质。数据范围是：$1 \leqslant A \leqslant B \leqslant 10^{15}$，$1 \leqslant N \leqslant 10^5$。

乍一看毫无头绪，仿佛怎么做都会超时……其实用容斥的想法就很容易了～

题目描述

HDU 4135 Co-prime

Given a number N, you are asked to count the number of integers between A and B inclusive which are relatively prime to N.
Two integers are said to be co-prime or relatively prime if they have no common positive divisors other than 1 or, equivalently, if their greatest common divisor is 1. The number 1 is relatively prime to every integer.

分析

先简化问题，可以先求出 $[1,b]$ 中与 n 互质个数，再减去 $[1,a-1]$ 中与 n 互质的个数。

要求出闭区间 $[1,b]$ 中有多少个数字与 $n$ 互质，我们可以考虑：求闭区间 $[1,b]$ 中有多少个数字与 $n$ 不互质。
与 $n$ 不互质的数字都是与 $n$ 有公共因子的数字。我们可以先将 $n$ 分解质因数，它的质因子肯定不会很多（似乎最多是12个？好像有个公式的……）。
首先考虑与 $n$ 有一个公共质因子的数字，可以方便地求出；但是这样有些数就会被重复考虑。具体地说，与 $n$ 有两个公共质因子的数字被计算了两次。所以我们要减去有两个公因子的数的个数。那么这样又会少算有三个公共质因子的数的数量……

很明显是容斥了，对于奇数个质因子的乘积的倍数的个数，将总数加上这个数量；如果是偶数个，就减去这个数量。
接下来用类似于状态压缩 DP 的方法枚举所有质因子的组合，判断奇偶、累加就可以了。
对于选择 $k$ 个质因子，那么我们就需要知道 $[1,b]$ 中 $lcm(P_i)$ 的倍数的个数。因为我们搞出来的 $p_i$ 都是质数，所以 $\displaystyle lcm(P_i) =\prod_{i=1}^{k} P_i$。直接乘起来就可以了～

代码

#define CLEAR(x) memset(x,0,sizeof(x))
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
const int maxn=50;
int T;
long long a,b,n;
vector<long long> yinzi;
inline long long llread(){
    long long ret=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9') {if (ch=='-') f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9') ret=(long long)ret*10+ch-'0',ch=getchar();
    return (long long)ret*f;
}
inline void MakeYinzi(long long n){
    yinzi.clear();
    for (long long i=2;i*i<=n;i++) if (n%i==0){
        yinzi.push_back(i);
        while (n%i==0) n/=i;
    }
    if (n>1) yinzi.push_back(n);
}
inline long long GetAnswer(long long x,long long n){
    long long ret=0;
    MakeYinzi(n);
    long long s=1<<yinzi.size();
    for (int i=1;i<s;i++){
        long long now=1,cnt=0;
        for (long long j=0;j<yinzi.size();j++) if (i&(1<<j)){
            now*=yinzi[j];
            cnt++;
        }
        if (cnt&1) ret+=x/now; else ret-=x/now;
    }
    return x-ret;
}
int main(){
    T=llread();
    for (int t=1;t<=T;t++){
        a=llread();b=llread();n=llread();
        long long ans=GetAnswer(b,n)-GetAnswer(a-1,n);
        printf("Case #%d: %lld\n",t,ans);
    }
    return 0;
}

拓展：ZOJ 3547

On Mars, there is a huge company called ACM (A huge Company on Mars), and it’s owned by a younger boss.

Due to no moons around Mars, the employees can only get the salaries per-year. There are n employees in ACM, and it’s time for them to get salaries from their boss. All employees are numbered from 1 to n. With the unknown reasons, if the employee’s work number is k, he can get k^4 Mars dollars this year. So the employees working for the ACM are very rich.

Because the number of employees is so large that the boss of ACM must distribute too much money, he wants to fire the people whose work number is co-prime with n next year. Now the boss wants to know how much he will save after the dismissal.

大概意思是：你要选出 $[1,n]$ 中与 $n$ 互质数字并且输出 $\displaystyle \sum_{i=1}^{r} i^4$。

几乎就是上面那题一样的题目，唯一的不同是需要用到四次方求和公式：

$$\displaystyle \sum_{i=1}^{n} i^4 = \frac {n(n+1)(2n+1)(3n^2+3n+1)} {30}$$

注意几个细节就可以了～代码：

#define CLEAR(x) memset(x,0,sizeof(x))
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
const int maxn=100;
const long long tt=1e9+7;
const int inv30=233333335;
int T;
long long n;
vector<long long> yinzi;
inline long long llread(){
    long long ret=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9') {if (ch=='-') f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9') ret=(long long)ret*10+ch-'0',ch=getchar();
    return (long long)ret*f;
}
inline void MakeYinzi(int n){
    yinzi.clear();
    for (int i=2;i*i<=n;i++) if (n%i==0){
        yinzi.push_back(i);
        while (n%i==0) n/=i;
    }
    if (n>1) yinzi.push_back(n);
}
inline long long Get4(long long x){
    long long ret= (long long)x*(x+1)%tt*(2*x+1)%tt*(3*x*x%tt+3*x%tt-1+tt)%tt;
    return ret;
}
inline long long Make(long long x){
    long long ret=Get4(n/x);
    return x*x%tt*x%tt*x%tt*ret%tt;
}
inline long long GetAnswer(long long n){
    MakeYinzi(n);
    long long ret=Get4(n);
    int s=1<<yinzi.size();
    for (int i=1;i<s;i++){
        long long now=1,cnt=0;
        for (int j=0;j<yinzi.size();j++) if (i&(1<<j)){
            now*=yinzi[j];
            cnt++;
        }
        if ((cnt&1)==0) ret=(ret+Make(now))%tt; else ret=(ret-Make(now)+tt)%tt;
    }
    return (ret*inv30)%tt;
}
int main(){
    T=llread();
    for (int t=1;t<=T;t++){
        n=llread();
        long long ans=GetAnswer(n);
        printf("%lld\n",ans);
    }
    return 0;
}