欧拉函数相关数论

欧拉函数

欧拉函数ϕ(n) 的百度百科定义是：“对正整数n，欧拉函数是小于n 的正整数中与n 互质的数的数目(ϕ(1)=1)。”不如将小于改成小于等于，倒是省去ϕ(1) 的特例。关于欧拉函数的证明，不多赘述（但仍强烈建议先看懂“证明”链接中的文章）。本文侧重于欧拉函数在编程竞赛中的应用，这里只将其中三个知识点提出来（其他公式皆可由这三个公式推导得出）：

1.ϕ(1)=1

2.ϕ(n)=n∏ir(1−1pi)其中n=∏1rpkii且pi为质数

3. ϕ(ab)=ϕ(a)ϕ(b)(积性函数性质，其中a、b互质)

素数打表

素数表较常与欧拉函数一同出现，除了一些特殊的数字外，若想要较快地求出某个数字的欧拉函数值，也需要先打出一张质数表来。打出小于n 的所有质数表的基本思想如下：

1. 用visit 数组记录值i 是否为合数，若为合数，则visit[i] = true，否则为false。用prime 数组储存素数的值，即所求素数表

2. 将i 从2 到n 遍历，如果当前值为素数（visit 值未被设置为true），则记录该数字到素数表中

3. 对于遍历的任何数字i，都用当前素数表中所有的素数与i 相乘，将乘积的visit 值设为true，表示该值为合数

4. 优化：当i 为当前素数表中某个素数的倍数时，跳出第三步的循环。这里需要说明一下，首先素数表中所有的数字为从小到大储存，且i 值也是从小到大遍历，举个例子：当我们遍历到i=15,prime=5(i%prime=0) 时，若不break，则会将visit[i*prime]=visit[75] 设为true，而往后遍历到i=25,prime=3 时，会再次把visit[75] 设为true，于是就出现了重复操作，这就是第四步所优化的地方。

第4 步优化的简单证明：

若i%primej=0，则对于primej+1，有：

i×primej+1=k×primej×primej+1=i′×primej

由primej+1>primej 可知：i′>i，由遍历次序可知，i′ 将随着大循环的i++ 在后面被遍历到，所以若出现i%primej=0 的情况，就不必再往后遍历primej+1 了，若觉得证明不太容易理解，可以先看代码再对照代码进行理解。

代码：

const int maxn = 10000 + 100;
int cnt;
int prime[maxn];
bool visit[maxn];
// n 表示要打的素数表最大值，cnt 表示素数的总数
void Prime(int n) {
cnt = 0;
memset(visit, 0, sizeof(visit));
for(int i = 2; i < n; ++i) {
if(!visit[i]) prime[cnt++] = i;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
visit[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}

n√ 求欧拉函数

要求欧拉函数，根据公式，就需要求出n 的所有素数因子，庆幸我们有这样一条性质：大于n√ 的素数因子最多只有1 个，因此我们只需要从2 找到n√，边找边除，就能得到最后一个素数。

证明：如果存在两个大于n√ 的素数因子，设分别为a 和b，则有a>n√,b>n√,n=abc>nc，其中c≥1，故矛盾。

由于可以在n√ 时间复杂度内求得n 的欧拉函数，根据一般题目的内存限制，我们可以打的素数表大小大概在106 左右，因此用这种方法，我们可以求得n≈1012 大小的欧拉函数。

代码

int euler(int n) {
if(n == 1) return 1;
int ret = n;
for(int j = 0; j < cnt && prime[j] * prime[j] <= n; ++j) {
if(n % prime[j] == 0) {
ret = ret / prime[j] * (prime[j] - 1);
while(n % prime[j] == 0) n /= prime[j];
}
}
if(n > 1) ret = ret / n * (n - 1);
return ret;
}

欧拉函数打表

除了打完素数表后根据公式求任意值的欧拉函数，其实也可在给素数打表的同时求出数字i 的欧拉函数值，相对于上面一种方法的限制，即不能求太大数字的欧拉函数，只能在内存限制范围内打表。

欧拉函数打表主要根据以下几点：

1. ϕ(1)=1

2. 若n 为素数，则ϕ(n)=n−1，显然对于任何素数n，所有小于它的正整数都与其互质。素数只有自己一个素数因子，代入上面的欧拉函数也可计算得n−1

3. 若i%prime=0 且n=i×prime，则

ϕ(n)=ϕ(k×prime2)=ϕ(k)ϕ(prime2)=ϕ(k)×prime2(1−1prime)=ϕ(k)ϕ(prime)×prime=ϕ(k×prime)×prime=ϕ(i)×prime

4. 若n=i×prime，则

ϕ(n)=ϕ(i×prime)=ϕ(i)ϕ(prime)=ϕ(i)×(prime−1)

由以上几点，就可以在打素数表的同时，打欧拉函数表了。

代码

const int maxn = 10000 + 100;
int cnt;
int prime[maxn], phi[maxn];
bool visit[maxn];
// n 表示要打的素数表最大值，cnt 表示素数的总数
void Prime(int n) {
cnt = 0;
phi[1] = 1;
memset(visit, 0, sizeof(visit));
for(int i = 2; i < n; ++i) {
if(!visit[i]) prime[cnt++] = i, phi[i] = i - 1;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
int k = i * prime[j];
visit[k] = true;
if(i % prime[j] == 0) {phi[k] = phi[i] * prime[j]; break;}
else phi[k] = phi[i] * (prime[j] - 1);
}
}
}

欧拉函数习题

基础知识部分完毕，下面是一些相关习题。

The Euler function

题目链接

HDU 2824: The Euler function

题意

多组数据，输入a,b，输出∑bi=aϕ(i) 的值，其中(2<a<b<3000000)。

题解

打表暴力加，注意一下用long long，前缀和预处理会爆空间

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long
#define For(i, a, b) for(int i = (a), _##i = (b); i <= _##i; ++i)
#define Rof(i, a, b) for(int i = (a), _##i = (b); i >= _##i; --i)
#define Mem(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(b))

const int maxn = 3000000 + 100;
int cnt, a, b;
LL ans;
int prime[maxn], phi[maxn];
bool visit[maxn];

void Prime(int n) {
cnt = 0;
phi[1] = 1;
Mem(visit, 0);
For(i, 2, n - 1) {
if(!visit[i]) prime[cnt++] = i, phi[i] = i - 1;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
int k = i * prime[j];
visit[k] = true;
if(i % prime[j] == 0) {phi[k] = phi[i] * prime[j]; break;}
else phi[k] = phi[i] * (prime[j] - 1);
}
}
}

int main() {
#ifdef LOCAL
freopen("test.txt", "r", stdin);
#endif // LOCAL
ios::sync_with_stdio(false);

Prime(maxn - 50);
while(scanf("%d%d", &a, &b) != EOF) {
ans = 0;
For(i, a, b) ans += phi[i];
printf("%I64d\n", ans);
}

return 0;
}

Happy 2006

题目链接

POJ 2773: Happy 2006

题意

多组数据，输入m,K，输出第K 小的与m 互质的数，其中(1≤m≤106),(1≤K≤108)。

题解

由gcd(m,n)=gcd(m,n+t×m)(t∈Z) 可得，所有与m 互质(gcd(m,n)=1) 的数以m 为周期变化，所以只需要求得这个周期(即ϕ(m))，再暴力在周期内搜一遍即可，要注意取膜等于0 的情况。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long
#define For(i, a, b) for(int i = (a), _##i = (b); i <= _##i; ++i)
#define Rof(i, a, b) for(int i = (a), _##i = (b); i >= _##i; --i)
#define Mem(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(b))

const int maxn = 3000000 + 100;
int cnt, m, k, ans, time;
int prime[maxn], phi[maxn];
bool visit[maxn];

void Prime(int n) {
cnt = 0;
phi[1] = 1;
Mem(visit, 0);
For(i, 2, n - 1) {
if(!visit[i]) prime[cnt++] = i, phi[i] = i - 1;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
int k = i * prime[j];
visit[k] = true;
if(i % prime[j] == 0) {phi[k] = phi[i] * prime[j]; break;}
else phi[k] = phi[i] * (prime[j] - 1);
}
}
}

int gcd(int a, int b) {
if(b == 0) return a;
return gcd(b, a % b);
}

int main() {
#ifdef LOCAL
freopen("test.txt", "r", stdin);
#endif // LOCAL
ios::sync_with_stdio(false);

Prime(maxn - 50);
while(scanf("%d%d", &m, &k) != EOF) {
time = k / phi[m];
k %= phi[m];
if(k == 0) k = phi[m],  --time;
For(i, 1, m) {
if(gcd(i, m) == 1) {
--k;
}
if(k == 0) {
ans = i;
break;
}
}
printf("%d\n", ans + time * m);
}

return 0;
}

Divisors

题目链接

POJ 2992: Divisors

题意

多组数据，输入n,k，输出Ckn 的约数个数，其中(0≤k≤n≤431), 且计算结果在long long 范围内。

题解

首先要求任意数字n 的约数个数，要先将n 分解质因数为n=∏ripkii，可以知道n 的所有约数都由这些数字组合相乘得到，任意质因数pi 都可以从0−ki 中任取一个次数去乘，因此n 的约数个数即∏ri(1+ki)。

其次，Ckn=∏ni=n−k+1ik!，所以可以将从1−431 所有数字的质因数字数表打出来，对所有质因数做前缀和（因为Ckn 的分子分母都由连续数字相乘），每次询问相减即可。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long
#define For(i, a, b) for(int i = (a), _##i = (b); i <= _##i; ++i)
#define Rof(i, a, b) for(int i = (a), _##i = (b); i >= _##i; --i)
#define Mem(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(b))

const int maxn = 450 + 100;
int cnt, n, k;
int prime[maxn];
LL sum[maxn][100], fac[maxn][100];
bool visit[maxn];

void Prime(int n) {
cnt = 0;
Mem(visit, 0);
For(i, 2, n - 1) {
if(!visit[i]) prime[cnt++] = i;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
visit[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}

void Cal_fac(int n) {
For(i, 1, n) {
int num = i;
For(j, 0, cnt - 1) {
while(num % prime[j] == 0) {
++fac[i][j];
num /= prime[j];
}
}
}
For(i, 1, n) {
For(j, 0, cnt - 1) {
sum[i][j] = sum[i - 1][j] + fac[i][j];
}
}
}

int main() {
#ifdef LOCAL
freopen("test.txt", "r", stdin);
#endif // LOCAL
ios::sync_with_stdio(false);

Prime(maxn - 50);
Cal_fac(maxn - 50);

while(scanf("%d%d", &n, &k) != EOF) {
LL ans = 1;
For(i, 0, cnt - 1) {
ans *= (1 + sum
[i] - sum[n - k][i] - sum[k][i]);
}
printf("%I64d\n", ans);
}

return 0;
}

Soldier and Number Game

题目链接

CodeForces 546D: Soldier and Number Game

题意

T 组数据，对于每一组数据的a、b，求a!b! 的素数因子个数。(1≤T≤106)(1≤b≤a≤5×106)

题解

对于106 组测试数据，肯定需要预处理，由于阶乘是连续数字相乘，所以可以预处理出从1 到n 中所有素数因子个数和，即n! 的素数因子个数，最后相减即可。

打表可得，在5×106 范围内的素数大概有34000 个素数，如果对1−5×106 内所有数字n 依次打表，计算量大概在15×1010，显然预处理会超时，需要对预处理进行优化。

这里只需要用到一点：如果n 能被素数prime 整除，则(n,n+prime) 之间所有数字都不能被prime 整除，这样便可通过前一个能够被prime 整除的数直接找到下一个能够被其整除的数，步长将随着prime 的增大而增大，而时间复杂度为n2+n3+n5+n7+⋯nn√=O(nlogn√)（最后一个素数若大于n√ 可直接求出）。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long
#define For(i, a, b) for(int i = (a), _##i = (b); i <= _##i; ++i)
#define Rof(i, a, b) for(int i = (a), _##i = (b); i >= _##i; --i)
#define Mem(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(b))

const int maxn = 5000000 + 100;
const int sq = sqrt((double)maxn);
int cnt;
int prime[maxn], phi[maxn];
LL num[maxn];
int tmp[maxn];
LL sum[maxn];
bool visit[maxn];

int read() {
char ch;
int ret = 0;
do {
ch = getchar();
} while(ch < '0' || ch > '9');
do {
ret = ret * 10 + ch - '0';
ch = getchar();
} while(ch >= '0' && ch <= '9');
return ret;
}

char str[20];
void write(LL n) {
if(n == 0) {putchar('0'); return ;}
int scnt = 0;
while(n != 0) {
str[++scnt] = n % 10 + '0';
n /= 10;
}
while(scnt != 0) putchar(str[scnt--]);
}

void Prime(int n) {
cnt = 0;
phi[1] = 1;
Mem(visit, 0);
For(i, 2, n - 1) {
if(!visit[i]) prime[cnt++] = i, phi[i] = i - 1;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
int k = i * prime[j];
visit[k] = true;
if(i % prime[j] == 0) {phi[k] = phi[i] * prime[j]; break;}
else phi[k] = phi[i] * (prime[j] - 1);
}
}
}

int main() {
#ifdef LOCAL
freopen("test.txt", "r", stdin);
#endif // LOCAL
ios::sync_with_stdio(false);

Prime(sq);
For(i, 1, maxn) tmp[i] = i;
For(i, 0, cnt - 1) {
for(int j = prime[i]; j <= maxn; j += prime[i]) {
while(tmp[j] % prime[i] == 0) {
++num[j];
tmp[j] /= prime[i];
}
}
}
For(j, 2, maxn - 1) {
if(tmp[j] != 1) ++num[j];
sum[j] = sum[j - 1] + num[j];
}
int t, a, b;
scanf("%d", &t);
while(t--) {
a = read(); b = read();
write(sum[a] - sum[b]);
putchar('\n');
}

return 0;
}

Longge’s problem

题目链接

POJ 2480: Longge’s problem

题意

多组输入，对于每一个输入数字N，求出∑Nigcd(i,N)，

题解

首先：

gcd(t,ab)=gcd(t,a)gcd(t,b)t∈Z∗,gcd(a,b)=1

由于ab 互质，所以a 和b 与任意正整数的最大公约数必然没有交集，上式成立，故gcd(t,n) 是积性函数，由积性函数性质可得：积性函数的和也是积性的，故∑nigcd(i,n) 也是积性函数。所以，若将n 分解为

n=∏irpkii(pi为质数)

则∑ingcd(i,n)=∏ir⎛⎝∑jpigcd(j,pkii)⎞⎠

其次，若gcd(a,b)=c，则gcd(ac,bc)=1，所以与a 的公约数值为c 的数字个数等于与ac 互质的数字个数相等，可以得到：

∑ingcd(i,n)=∑imfaciϕ(nfaci)faci为n的约数，m为n的约数个数

代入上式继续推导可得：

∑ingcd(i,n)=∏ir⎛⎝∑jmfacjϕ(pkiifacj)⎞⎠

由于素数的k 次幂的所有约数分别为该素数的0−k 次幂，所以

∑ingcd(i,n)=∏ir⎛⎝∑jkipjiϕ(pkiipji)⎞⎠=∏ir⎛⎝∑jkipjiϕ(pki−ji)⎞⎠

若推导到这里便开始写程序，容易超时，还需要一步计算：若n=primek，则ϕ(n)=primek(1−1prime)=primek−primek−1，所以：

∑ingcd(i,n)=∏ir⎛⎝∑jki−1pji(pki−ji−pki−j−1i)+pkii⎞⎠=∏ir⎛⎝∑jki−1(pkii−pki−1i)+pkii⎞⎠=∏ir(ki(pkii−pki−1i)+pkii)

注意这里i 从1 到r，而j 从0 到ki，最后当j=ki 时，pkiiϕ(1)，而ϕ(1) 的值为特判，不能由欧拉公式直接计算得到，故应单独分离出来进行计算。

过题代码

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <climits>
#include <cfloat>
#include <cstring>
#include <string>
#include <vector>
#include <list>
#include <queue>
#include <stack>
#include <map>
#include <set>
#include <algorithm>
using namespace std;

#define LL long long
#define For(i, a, b) for(int i = (a), _##i = (b); i <= _##i; ++i)
#define Rof(i, a, b) for(int i = (a), _##i = (b); i >= _##i; --i)
#define Mem(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(b))

const int maxn = 100000 + 100;
int cnt, fac;
LL ans, n;
int prime[maxn];
bool visit[maxn];

void Prime(int n) {
cnt = 0;
Mem(visit, 0);
For(i, 2, n - 1) {
if(!visit[i]) prime[cnt++] = i;
for(int j = 0; j < cnt && i * prime[j] < n; ++j) {
visit[i * prime[j]] = true;
if(i % prime[j] == 0) break;
}
}
}

int main() {
#ifdef LOCAL
freopen("test.txt", "r", stdin);
#endif // LOCAL
ios::sync_with_stdio(false);

Prime(maxn - 50);
while(scanf("%I64d", &n) != EOF) {
ans = 1;
int sq = sqrt((double)n);
for(int i = 0; prime[i] <= sq; ++i) {
fac = 0;
while(n % prime[i] == 0) {
++fac;
n /= prime[i];
}
ans *= fac * (pow(prime[i], fac) - pow(prime[i], fac - 1)) + pow(prime[i], fac);
}
if(n != 1) ans *= 2 * n - 1;
printf("%I64d\n", ans);
}

return 0;
}