欧拉函数相关数论

博文

---

欧拉函数

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

$$\phi(1)=1$$
2. $$\phi(n)=n\prod_i^r\left(1-\frac{1}{p_i}\right)\\其中\quad n=\prod_1^rp_i^{k_i}\quad 且p_i为质数$$
3. $$\phi(ab)=\phi(a)\phi(b)\quad(积性函数性质，其中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\%prime_j=0$，则对于 $prime_{j+1}$，有：

$$i\times prime_{j+1}=k\times prime_j\times prime_{j+1}=i'\times prime_j$$
由 $prime_{j+1}>prime_j$ 可知：$i'>i$，由遍历次序可知，$i'$ 将随着大循环的 i++ 在后面被遍历到，所以若出现 $i\%prime_j=0$ 的情况，就不必再往后遍历 $prime_{j+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;
        }
    }
}

$\sqrt n$ 求欧拉函数

要求欧拉函数，根据公式，就需要求出 $n$ 的所有素数因子，庆幸我们有这样一条性质：大于 $\sqrt n$ 的素数因子最多只有 $1$ 个，因此我们只需要从 $2$ 找到 $\sqrt n$，边找边除，就能得到最后一个素数。
证明：如果存在两个大于 $\sqrt n$ 的素数因子，设分别为 $a$ 和 $b$，则有 $a>\sqrt n,b>\sqrt n,n=abc>nc$，其中 $c\geq 1$，故矛盾。
由于可以在 $\sqrt n$ 时间复杂度内求得 $n$ 的欧拉函数，根据一般题目的内存限制，我们可以打的素数表大小大概在 $10^6$ 左右，因此用这种方法，我们可以求得 $n\approx 10^{12}$ 大小的欧拉函数。

代码

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. $\phi(1)=1$;
2. 若 $n$ 为素数，则 $\phi(n)=n-1$，显然对于任何素数 $n$，所有小于它的正整数都与其互质。素数只有自己一个素数因子，代入上面的欧拉函数也可计算得 $n-1$;
3. 若 $i\%prime=0$ 且 $n=i\times prime$，则

$$\phi(n)=\phi(k\times prime^2)=\phi(k)\phi(prime^2)=\phi(k)\times prime^2\left(1-\frac{1}{prime}\right)\\=\phi(k)\phi(prime)\times prime=\phi(k\times prime)\times prime=\phi(i)\times prime$$
4. 若 $n=i\times prime$，则
$$\phi(n)=\phi(i\times prime)=\phi(i)\phi(prime)=\phi(i)\times(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$，输出 $\sum_{i=a}^b\phi(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;
LL 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);
        }
    }
    For(i, 1, n - 1) {
        phi[i] += phi[i - 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) {
        printf("%I64d\n", phi[b] - phi[a - 1]);
    }
    return 0;
}

Happy 2006

题目链接

POJ 2773: Happy 2006

题意

多组数据，输入 $m,K$，输出第 $K$ 小的与 $m$ 互质的数，其中 $(1\leq m\leq10^6),(1\leq K\leq10^8)$。

题解

由 $gcd(m,n)=gcd(m,n+t\times m)(t\in Z)$ 可得，所有与 $m$ 互质 $(gcd(m,n)=1)$ 的数以 $m$ 为周期变化，所以只需要求得这个周期 $(即\phi(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$，输出 $C_n^k$ 的约数个数，其中 $(0\leq k\leq n\leq431)$, 且计算结果在 long long 范围内。

题解

首先要求任意数字 $n$ 的约数个数，要先将 $n$ 分解质因数为 $n=\prod_i^rp_i^{k_i}$，可以知道 $n$ 的所有约数都由这些数字组合相乘得到，任意质因数 $p_i$ 都可以从 $0-k_i$ 中任取一个次数去乘，因此 $n$ 的约数个数即 $\prod_i^r(1+k_i)$。
其次，$C_n^k=\frac{\prod_{i=n-k+1}^ni}{k!}$，所以可以将从 $1-431$ 所有数字的质因数字数表打出来，对所有质因数做前缀和（因为 $C_n^k$ 的分子分母都由连续数字相乘），每次询问相减即可。

过题代码

#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[n][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$，求 $\frac{a!}{b!}$ 的素数因子个数。$(1\leq T\leq10^6)(1\leq b\leq a\leq5\times10^6)$

题解

对于 $10^6$ 组测试数据，肯定需要预处理，由于阶乘是连续数字相乘，所以可以预处理出从 $1$ 到 $n$ 中所有素数因子个数和，即 $n!$ 的素数因子个数，最后相减即可。
打表可得，在 $5\times10^6$ 范围内的素数大概有 $34000$ 个素数，如果对 $1-5\times10^6$ 内所有数字 $n$ 依次打表，计算量大概在 $15\times10^{10}$，显然预处理会超时，需要对预处理进行优化。
这里只需要用到一点：如果 $n$ 能被素数 $prime$ 整除，则 $(n,n+prime)$ 之间所有数字都不能被 $prime$ 整除，这样便可通过前一个能够被 $prime$ 整除的数直接找到下一个能够被其整除的数，步长将随着 $prime$ 的增大而增大，而时间复杂度为 $\frac{n}{2}+\frac{n}{3}+\frac{n}{5}+\frac{n}{7}+\cdots\frac{n}{\sqrt n}=O(n\log \sqrt n)$（最后一个素数若大于$\sqrt 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$，求出 $\sum_i^Ngcd(i,N)$，

题解

首先：

$$gcd(t,ab)=gcd(t,a)gcd(t,b)\quad \quad t\in Z^*,gcd(a,b)=1$$
由于 $ab$ 互质，所以 $a$ 和 $b$ 与任意正整数的最大公约数必然没有交集，上式成立，故 $gcd(t,n)$ 是积性函数，由积性函数性质可得：积性函数的和也是积性的，故 $\sum_i^ngcd(i,n)$ 也是积性函数。所以，若将 $n$ 分解为
$$n=\prod_i^rp_i^{k_i}\quad\quad(p_i 为质数)$$
则
$$\sum_i^ngcd(i,n)=\prod_i^r\left(\sum_j^{p_i}gcd(j,p_i^{k_i})\right)$$
其次，若 $gcd(a,b)=c$，则 $gcd(\frac{a}{c},\frac{b}{c})=1$，所以与 $a$ 的公约数值为 $c$ 的数字个数等于与 $\frac{a}{c}$ 互质的数字个数相等，可以得到：
$$\sum_i^ngcd(i,n)=\sum_i^mfac_i\phi\left(\frac{n}{fac_i}\right)\quad\quad fac_i 为n 的约数，m为n的约数个数$$
代入上式继续推导可得：
$$\sum_i^ngcd(i,n)=\prod_i^r\left(\sum_j^mfac_j\phi\left(\frac{p_i^{k_i}}{fac_j}\right)\right)$$
由于素数的 $k$ 次幂的所有约数分别为该素数的 $0-k$ 次幂，所以
$$\sum_i^ngcd(i,n)=\prod_i^r\left(\sum_j^{k_i}p_i^j\phi\left(\frac{p_i^{k_i}}{p_i^j}\right)\right)=\prod_i^r\left(\sum_j^{k_i}p_i^j\phi\left(p^{k_i-j}_i\right)\right)$$
若推导到这里便开始写程序，容易超时，还需要一步计算：若 $n=prime^k$，则 $\phi(n)=prime^k(1-\frac{1}{prime})=prime^k-prime^{k-1}$，所以：
$$\sum_i^ngcd(i,n)=\prod_i^r\left(\sum_j^{k_i-1}p^j_i\left(p_i^{k_i-j}-p_i^{k_i-j-1}\right)+p_i^{k_i}\right)\\ =\prod_i^r\left(\sum_j^{k_i-1}\left(p_i^{k_i}-p_i^{k_i-1}\right)+p_i^{k_i}\right)=\prod_i^r\left(k_i\left(p_i^{k_i}-p_i^{k_i-1}\right)+p_i^{k_i}\right)$$
注意这里 $i$ 从 $1$ 到 $r$，而 $j$ 从 $0$ 到 $k_i$，最后当 $j=k_i$ 时，$p_i^{k_i}\phi(1)$，而 $\phi(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;
}