能力综合提升-数学-1

0x01 素数

线性筛素数

首先是线性筛吧！这个必须要背的。

注意：在线性筛中，每一个合数都是被最小的质因子筛掉。

时间复杂度 $O(n)$。

vector<int> prime;
bitset<Max> notprime;

void sift(const int Max)
{
    F(i, 2, Max)
    {
        if(!notprime[i]) prime.push_back(i);
        for(int j : prime)
        {
            if(i * j > Max) break;
            notprime[i * j] = true;
            if(i % j == 0) break; 
        }
    }
}

Miller-Rabin

判断一个数是不是素数。

1. 暴力试除法 $O(n),O(\sqrt{n})$，不可接受。
2. Fermat 小定理：$p$ 为素数时，$a^p\equiv a(\mathrm{mod}p)$ —— 然而逆命题并不一定成立，这些素数被称为 Carmichael 数（$561=3\times 11\times 17$）。在 ${10}^9$ 以内，这样的数有 646 个，因此不能满足需求。
3. Miller-Rabin 素性测试。

以下内容来自 warzone 的 Blog。

Miller-Rabin 素性测试是对 Fermat 小定理的优化。具体的，我们有二次探测定理：

若 $p$ 为奇素数，则 $x^2\equiv 1(\mathrm{mod}p)$ 的解为 $x\equiv \pm 1(\mathrm{mod}p)$

PPT 上是这么说的：如果 $n$ 为素数，取 $a<n$，设 $n-1=d\times 2^r$，则要么 $a^d\equiv 1(\mathrm{mod}n)$，要么 $\mathrm{\exists }0\le i<r,\mathrm{s}.\mathrm{t}.a^{d\times 2^i}\equiv -1(\mathrm{mod}n)$。其证明见 Wikipedia。

具体的，若 $p$ 是奇数，则 $p-1$ 是偶数：

1. 选取一个底数 $a$，初始化指数 $d=p-1$，判断 $a^{d}modp$ 是否为 $1$，若否则非素数。
2. 重复 $d\leftarrow {\displaystyle \frac{d}{2}}$ 并判断，直到 $d$ 为奇数。
3. 最后，若 $a^d\not\equiv \pm 1(\mathrm{mod}p)$，则 $p$ 非素数，否则 $p$ 大概率是一个素数。

bool check(const int a, const ll p)
{
    ll d = p - 1, res = fastpow(a, d, p);
    if(res != 1) return false;
    while(!(d & 1))
    {
        res = fastpow(a, d >>= 1, p);
        if(res == p - 1) return true;
        else if(res != 1) return false;
    }
    return true;
}

初学者写 Miller-Rabin 时可能会犯的错误：

1. 没有特判 $2^{p-1}\equiv -1(\mathrm{mod}p)$ 的情况。
2. 先计算 $a^{d}modp$，再检查 $d$ 的奇偶性。

若随机选取 $k$ 个底数，Miller-Rabin 的错误率为 ${\displaystyle \frac{1}{4^k}}$。也就是说，选取 $20\sim 40$ 个底数，正确率是可以接受的。

对于竞赛的特殊环境，有

1. 判断 unsigned int 范围内，选取 2, 7, 61 三个底数即可。
2. 判断 unsigned long long 范围内，选取 2, 325, 9375, 28178, 450775, 9780504, 1795265022 七个底数即可。
3. 判断 unsigned long long 范围内，选取前十二个素数（到 37 为止）作为底数即可，OEIS。

注：不是哥们，你这个写的太他妈丑陋了，还是换一个抄吧。

注：fastpow 要注意别炸 long long，用 __int128。

inline ll fastpow(ll a, ll b, const ll mod)
{
    ll ans = 1;
    while (b)
    {
        if (b & 1)
            ans = (__int128)ans * a % mod;
        a = (__int128)a * a % mod;
        b >>= 1;
    }
    return ans;
}
bool check(const ll a, const ll n)
{
    ll d = n - 1;
    while(d)
    {
        ll res = fastpow(a, d, n);
        if(res != 1 && res != n - 1)
            return false;
        if((d & 1) || res == n - 1)
            return true;
        d >>= 1;
    }
    return true;
}
vector<int> MRTest = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
bool Miller_Rabin(const ll p)
{
    if(p > 40)
    {
        for(int i : MRTest)
            if(!check(i, p)) return false;
        return true;
    }
    else
    {
        for(int i : MRTest)
            if(i == p) return true;
        return false;
    }
}

Pollard-Rho

质因数分解。

1. 暴力枚举依然是 $O(\sqrt{n})$ 的。

2. 利用生日悖论减少枚举量。

生日悖论：从 $n$ 个数中随机选取 $k$ 个两两之差全部不等于 $c$ 的概率约为 $O({\displaystyle \frac{1}{n^{k^2}}})$。

所以：找到 $k$ 个数，检查这 $k$ 个数两两的差与 $n$ 的 $gcd$。

伪随机函数 $f_i=f_{i-1}^2+c$

它具有的优秀的性质：如果 $gcd(f_i-f_j,n)>1$，则 $gcd(f_{i+1}-f_{j+1},n)=gcd((f_i-f_j)(f_i+f_j),n)>1$。

mt19937_64 rand64(time(nullptr));
ll rnd(const ll l, const ll r) // 左闭右开
{
    return l + (rand64() % (r - l));
}
ll f(const ll x, const ll c, const ll p)
{
    return (x * x + c) % p;
}
ll Pollard_Rho(const ll n)
{
    ll x = rnd(0, n), y = x;
    // 关于这里 y 和下面 k 的含义见下面正文
    ll c = rnd(3, n), d = 1;
    // 当 c = 0 或 c = 2 的情况应当避免
    // 来自算法导论，书上也没有具体解释原因
    for(int i = 1, k = 2; d == 1; i ++)
    {
        x = f(x, c, n);
        d = gcd(qabs(y - x), n);
        if(i == k)
        {
            y = x;
            k <<= 1;
        }
    }
    return d;
}

但是有一个问题：这个伪随机数的函数会包含死循环————循环节的产生不难理解：在模 $n$ 意义下，整个函数的值域只能是 $0,1,2,\cdots ,n-1$。总有一次，函数在迭代的过程中会带入之前的值，得到相同的结果。生成数列的轨迹很像一个希腊字母 $\rho$ ，所以整个算法的名字叫做 Pollard-Rho。

这个问题可以基于 Floyd 算法优化。具体的，我们希望判断环的存在，通俗来讲就是“龟兔赛跑”。假设环的长度为 $c$，在环内有 $f_i=f_{i+kc}$。乌龟每次迭代一次，兔子每次迭代两次，如果二者相同，则必有 $f_{2i}=f_{i+kc}$，这样就找到了环。

ll Pollard_Rho(const ll n)
{
    if(n == 4) // 初始跳两步没法分解 4，需要特判
        return 2;
    ll x = rnd(0, n), y = x;
    ll c = rnd(3, n), d = 1;
    x = f(x, c, n);
    y = f(f(y, c, n), c, n);
    while(x != y)
    {
        x = f(x, c, n);
        y = f(f(y, c, n), c, n);
        d = gcd(abs(x - y), n);
        if(d > 1)
            return d;
    }
    return n;
}

考虑做 $gcd$ 的时间复杂度是 $O(\log n)$ 的，如果多次使用时间复杂度会大到爆炸。其实，每次计算 $gcd$ 是不必要的，可以计算一段连续的乘积，然后再做 $gcd$ 也是一样的效果（倍增的思想）。

1. 枚举 $i=2^j,2^j+1,\cdots ,2^{j+1}-1$，计算 $gcd(\prod (f_i-f_{2^j}),p)$。
2. 但是这样长度太长才检查一次，所以加一个每 $t$ 次乘积之后就检查一次的条件。一般 $t=128$（此处无任何道理说是）。

注：做乘法要防止炸 long long，必须使用 __int128。

ll mul(ll a, ll b, ll p)
{
    return __int128(a) * b % p;
}
mt19937_64 rand64(time(nullptr));
ll rnd(const ll l, const ll r)
{
    return l + (rand64() % (r - l));
}
ll f(const ll x, const ll c, const ll p)
{
    return (mul(x, x, p) + c) % p;
}
int cnt;
ll Pollard_Rho(ll n)
{
    if(n == 4)
        return 2;
    ll x = rnd(0, n), y = x;
    ll c = rnd(3, n), d = 1;
    x = f(x, c, n);
    y = f(f(y, c, n), c, n);
    for(int lim = 1; x != y; lim = qmin(128, lim << 1))
    {
        ll cnt = 1;
        F(i, 1, lim)
        {
            ll tmp = mul(cnt, qabs(x - y), n);
            if(tmp == 0)
                break;
            cnt = tmp;
            x = f(x, c, n);
            y = f(f(y, c, n), c, n);
        }
        d = gcd(cnt, n);
        if(d > 1)
            return d;
    }
    return n;
}   

---

0x02 最大公约数

辗转相除法

首先是最简单的欧几里得算法！

ll gcd(ll a, ll b)
{
    return b == 0 ? a : gcd(b, a % b);
}

不过时间复杂度是 $O(\log V)$。

当我们需要访问 gcd() 很多次时，就需要用到其他的算法了。

一种基于值域预处理的快速 gcd

以下内容来自 moongazer 的 Blog。

这是一种 $O(V)$ 时间预处理 $O(1)$ 询问求 gcd 的算法。

记：$V$ 为值域，$\mathbb{P}$ 为素数集，集合 $a_1,a_2,\cdots ,a_n$ 是 $x$ 的分解 $⟺x=\prod _{i=1}^n{a}_i$。

定理一：可以把值域内的每个 $x$ 分解成 $a,b,c$ 且必有 $a,b,c$ 要么 $\le \sqrt{x}$ 要么 $\in \mathbb{P}$（称为合法分解）。

反证。不失一般性，令 $a\le b\le c$。若 $c\notin \mathbb{P}\wedge c>\sqrt{x}$，则必存在 $c$ 的分解 $d,e$ 且有 $d\le e$ 和 $d\le \sqrt{x}$。

而 $a\times b=\frac{x}{c}<\sqrt{x}$，于是也有 $x$ 的分解 $d,a\times b,e$。若 $e>\sqrt{x}$ 则重复上述过程直到满足条件。

定理二：对于 gcd(x, y)，考虑 $x$ 的合法分解 $a,b,c$，分别考虑它们对答案的贡献。$a$ 的贡献为 $v_a=gcd(a,y)$，$b$ 的贡献为 $v_b=gcd(b,\frac{y}{v_a})$ （因为 $v_a$ 这个因子已经被计算过，不能重复计算），$c$ 的贡献同理，$v_c=gcd(c,\frac{y}{v_a{v}_b})$。

我们发现实现的难点在于如何在 $O(V)$ 时间内进行分解，询问部分较为容易。

记最小质因子函数 $f(x)$ 为 $x$ 的最小质因子。

如何分解？对于 $x\ge 2$，找到 $x$ 的最小质因子 $p=f(x)$ 和 $\frac{x}{p}$ 的一个合法分解 $a_0,b_0,c_0$，则 $x$ 的一个合法分解为 $\mathrm{sort}{a}_0\times p,b_0,c_0$。

试着证明它：

1. $x\in \mathbb{P}$ 时显然成立，分解为 $1,1,x$；
2. $p\le \sqrt[4]{x}$ 时有 $a_0\le \sqrt[3]{\frac{x}{p}}$ 则 $a_0\times p\le x^{\frac{1}{3}}{p}^{\frac{2}{3}}\le x^{\frac{1}{3}}{x}^{\frac{1}{6}}=\sqrt{x}$；
3. $p>\sqrt[4]{x}$ 时
   1. $a_0=1$，有 $a_0\times p=p\le \sqrt{x}$；
   2. $a_0\ne 1$，则 $x\notin \mathbb{P}$，$q=f(\frac{x}{p})\ge f(x)=p$，则 $p\le q\le a_0\le b_0\le c_0$，从而 $p\times a_0\times b_0\times c_0>(\sqrt[4]{x}{)}^4=x$ 与 $p\times a_0\times b_0\times c_0=x$ 矛盾，故此情况不存在。

constexpr int maxm = 1e6 + 10;
constexpr int maxk = 1e3 + 10;
struct Fac
{
    int a, b, c;
    Fac() : a(1), b(1), c(1) {}
    int operator[](int idx) // 方便使用 for 来调用每个元素
    {
        switch(idx)
        {
            case 1: return a;
            case 2: return b;
            case 3: return c;
            default: break;
        }
        return 0;
    }
    void rec(Fac t, int p) // 赋值并让 a *= p，并排序
    {
        a = t.a * p, b = t.b, c = t.c;
        sort();
    }
    void sort() // 注意到 bc 一定是顺序，只需要看 a 的位置即可
    {
        if(a > b) qswap(a, b);
        if(b > c) qswap(b, c);
    }
} fac[maxm]; //用于存放分解
bitset<maxm> isPrime; // bitset 存是否是质数
int Pri[maxm], cnt; // 质数集 P
int pre[maxk][maxk]; // 预处理的 sqrt(V) 内的 gcd
const int V = 1e6, T = 1e3; // T = sqrt(V)
void Pre() // 线性筛预处理
{
    isPrime.flip(); // 全部标为质数
    F(i, 2, V) // 线性筛
    {
        if(isPrime[i])
        {
            fac[i].c = i;
            Pri[++ cnt] = i;
        }
        for(int j = 1, k = Pri[j] * i; k <= V; j ++, k = Pri[j] * i)
        {
            isPrime[k] = false;
            fac[k].rec(fac[i], Pri[j]); // 把质数的分解偷过来，并且 Pre[j] 一定是最小质因子 p
            if(i % Pri[j] == 0)         // 最小的保证是因为 i % Pri[j] == 0 时就 break 掉了
                break;
        }
    }
    // 预处理 gcd 数组
    // 因为由辗转相除法，a > b 时有 gcd(a, b) = gcd(b, a % b)，这可以预处理实现
    F(i, 0, sqrt(V))
        pre[0][i] = pre[i][0] = i;
    F(i, 1, sqrt(V))
        F(j, 1, i)
            pre[i][j] = pre[j][i] = pre[j][i % j];
}
int gcd(int a, int b)
{
    int ans = 1;
    F(i, 1, 3)
    {
        int tmp;
        if(fac[a][i] > sqrt(V)) // 查看 a 的分解的第 i 个和 sqrt(V) 的关系
        {
            if(b % fac[a][i])
                tmp = 1; // 互质，没有贡献
            else
                tmp = fac[a][i]; // 整除，贡献就是本身
        }
        else
            tmp = pre[fac[a][i]][b % fac[a][i]]; // 否则两边都是 sqrt(V) 内的内容可以预处理解决
        b /= tmp; // 防止算重
        ans *= tmp; // 计算贡献
    }
    return ans;
}

例题

求出 $\mathrm{\forall }k=1,2,\cdots ,n$，在给定的 $a_i$ 中挑出 $k$ 个数，最大的 $gcd$。

$n\le {10}^4,a_i\le {10}^6$。

一个相当经典的结论是，直接去求 $gcd$ 是没有必要的，枚举因子反而是更简单也更正确的。

对于每个 $a_i$，$O(\sqrt{n})$ 暴力枚举因子把它们存到一个桶里，并且显然有 $an{s}_k\ge an{s}_{k+1}$。

那你从上往下扫一遍就好啦。

int main()
{
    int n = readint();
    F(i, 1, n)
    {
        read(a);
        F(j, 1, sqrt(a))
            if(a % j == 0)
                b[j] ++, b[a / j] += (a != j * j);
    }
    for(int i = 1, j = 1e6; i <= n; i ++)
    {
        while(b[j] < i)
            j --;
        write('\n', j);
    }
    return 0;
}

例题

$n$ 次询问，求出所有可能的 $x$ 的个数，使 $gcd(x,a_0)=a_1\wedge \mathrm{lcm}(x,b_0)=b_1$。

$n\le 2000,1\le a_0,a_1,b_0,b_1\le 2\times {10}^9$。

你看到这个东西很自然的就要想到

$$\{\begin{array}{l}gcd(x,a_0)=a_1⟹(xmod{a}_1=a_{0}mod{a}_1=0)\wedge (gcd(\frac{x}{a_1},\frac{a_0}{a_1})=1)\\ \mathrm{lcm}(x,b_0)=b_1⟹(b_{1}modx=b_{1}mod{b}_0=0)\wedge (gcd(\frac{b_1}{x},\frac{b_1}{b_0})=1)\end{array}$$

那你直接 $O(\sqrt{V})$ 枚举 $x$ 然后依次检查这四个条件即可。总时间复杂度 $O(n\sqrt{V})$ 在 ${10}^9$ 级别但是跑的飞快。

bool check(int c2, int d2, int a1, int b1, int x)
{
    int d1 = b1 / x;
    if(x % a1)
        return false;
    int c1 = x / a1;
    if(gcd(c1, c2) > 1)
        return false;
    if(gcd(d1, d2) > 1)
        return false;
    return true;
}
int workwork()
{
    int a0, a1, b0, b1;
    read(a0, a1, b0, b1);
    if(a0 % a1)
        return 0;
    int c2 = a0 / a1;
    if(b1 % b0)
        return 0;
    int d2 = b1 / b0;
    int ans = 0;
    F(i, 1, sqrt(b1))
    {
        if(b1 % i)
            continue;
        ans += check(c2, d2, a1, b1, i);
        if(i * i < b1)
            ans += check(c2, d2, a1, b1, b1 / i);
    }
    return ans;
}

---

0x03 欧拉函数

$\phi (n)=\sum _{i=1}^n[gcd(n,i)=1]$。

性质

• $n\in \mathbb{P}⟹\phi (n)=n-1$

• 是积性函数。即 $\mathrm{\forall }a,b,gcd(a,b)=1⟹\phi (ab)=\phi (a)\phi (b)$

• $n$ 为奇数时 $\phi (2n)=\phi (n)$

• $n=\sum _{d\mid n}\phi (d)$

• $n=p^k,p\in \mathbb{P},\phi (n)=p^k-p^{k-1}$

• $\phi (n)=n\prod _{i=1}^s{\displaystyle \frac{p_i-1}{p_i}}$，$p_i$ 是质因数

• $\phi (gcd(m,n))\phi (mn)=gcd(m,n)\phi (m)\phi (n)$

• $n>1⟹2\mid \phi (n)$

部分证明见 OI-Wiki。

单点求 $\phi (n)$

主要利用性质 $\phi (n)=n\prod _{i=1}^s{\displaystyle \frac{p_i-1}{p_i}}$（$p_i$ 是质因数）进行 $O(\sqrt{n})$ 的质因数分解。

可以用 Pollard-Rho 优化。

int Phi(int n)
{
    int ans = n;
    for(int i = 2; i * i <= n; i ++)
        if(n % i == 0) // 对于每个质因数
        {
            ans = ans / i * (i - 1);
            while(n % i == 0) // 除掉所有的当前质因数
                n /= i;
        }
    if(n > 1) // 分解最后剩下单个质数
        ans = ans / n * (n - 1);
    return ans;
}

线性筛 $\phi (n)$

只需要在线性筛素数的基础上改进即可。

还记得线性筛的原理吗？每一个合数都是被最小的质因子筛掉。

记 $p_1$ 为 $n$ 的最小质因子，有 $n=p_1\times n^{\prime }$。

分两个情况讨论：

• $n^{\prime }mod{p}_1=0$，那么 $n^{\prime }$ 包含了 $n$ 的所有质因子。

$\begin{array}{rl}\phi (n)& =n\times \prod {\displaystyle \frac{p_i-1}{p_i}}\\ & =p_1\times n^{\prime }\times \prod {\displaystyle \frac{p_i-1}{p_i}}\\ & =p_1\times \phi (n^{\prime })\end{array}$

• $n^{\prime }mod{p}_1\ne 0$，那么 $gcd(p_1,n^{\prime })=1$，由欧拉函数的积性有

$\begin{array}{rl}\phi (n)& =\phi (p_1\times n^{\prime })\\ & =\phi (p_1)\times \phi (n^{\prime })\\ & =(p_1-1)\times \phi (n^{\prime })\end{array}$

其中 $n^{\prime }<n$ 即 $\phi (n^{\prime })$ 已经被算出。

vector<int> prime;
bitset<Max> notprime;
int Phi[Max];

void sift(const int Max)
{
    Phi[1] = 1; // 特判
    F(i, 2, Max)
    {
        if(!notprime[i]) 
        {
            prime.push_back(i);
            Phi[i] = i - 1; // 质数的情况
        }
        for(int j : prime)
        {
            if(i * j > Max) break;
            notprime[i * j] = true;
            if(i % j == 0)
            {
                Phi[i * j] = Phi[i] * j; // 情况 1，直接乘 p_1
                break;
            }
            Phi[i * j] = Phi[i] * Phi[j]; // 情况 2，乘 Phi(p_1) 也就是 (p_1 - 1)
        }
    }
}

例题

$n\times n$ 的方阵，求视线能看到的人的个数。

忽略几个特殊的点，我们注意到，可以被看到的人必须满足 $gcd(x,y)=1$。 证明是显然的，若 $d=gcd(x,y)\ne 1$，则 $({\displaystyle \frac{x}{d}},{\displaystyle \frac{y}{d}})$ 必然能把 $(x,y)$ 挡住，所以我们只需要求出范围内互质的对数即可。

由于把 $(1,1)$ 看成坐标原点，因此 $x,y\le n-1$ 而非 $n$。

注意到 $i,j$ 的对称性，只需要计算 $i>j$ 的情况。剩下的 $(0,1),(1,0),(1,1)$ 最后补上即可。

$\begin{array}{rl}ans& =2+\sum _{i=1}^{n-1}\sum _{j=1}^{n-1}[gcd(i,j)=1]\\ & =2+\sum _{i=1}^{n-1}[gcd(i,i)=1]+2\times \sum _{i=1}^{n-1}\sum _{j=1}^{i-1}[gcd(i,j)=1]\\ & =3+2\times \sum _{i=1}^{n-1}\phi (i)\end{array}$

求 $\phi (n)$ 前 $n-1$ 项和即可。时间复杂度 $O(n)$。

例题

求 $\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)\in \mathbb{P}]$。

$\begin{array}{rl}ans& =\sum _{p\in \mathbb{P}}\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=p]\\ & =\sum _{p\in \mathbb{P}}\sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{p}\rfloor }[gcd(i,j)=1]\\ & =\sum _{p\in \mathbb{P}}(-1+2\times \sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\sum _{j=1}^i[gcd(i,j)=1])\\ & =\sum _{p\in \mathbb{P}}(-1+2\times \sum _{i=1}^{\lfloor \frac{n}{p}\rfloor }\phi (i))\end{array}$

$O(n)$ 预处理 $\phi (n)$ 及 $\sum _{i=1}^n\phi (i)$，枚举每个质数即可。时间复杂度 $O(n)$。

例题

求 $\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)$。

$\begin{array}{rl}ans& =\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^n[gcd(i,j)=k]\times k\\ & =\sum _{k=1}^{n}k\times \sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=k]\end{array}$

对于 $\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=k]$ 的处理，套路同上。时间复杂度 $O(n)$。