The secret of being a bore is to tell everything

0%

质因数分解

质因数分解

接着上一章的素数判定,这一章我们要讨论质因数分解的算法。

算数基本定理,又称唯一分解定理:每个大于$1$的自然数,都可以写成质数次方的积,且这些质数从小到大排列后,仅存在这一种分解方式。

换句话说,每个$n>1$,都有唯一一种分解方式$\prod_{k} p_k^{\alpha_k} = n$。

朴素方法

对于所有$p|n$,我们不断试除就可以了,与质数判定同理,我们不需要试超过$\sqrt{n}$的质数,因为这样的质数一定就是$n$本身,于是有代码:

1
map<ll, ll> factorize(ll x) {
2
    map<ll, ll> pcnt;
3
    for (ll i = 2; i * i <= x; i++) {
4
        if (x % i == 0) {
5
            int cnt = 0;
6
            while (x % i == 0) cnt++, x /= i;
7
            pcnt[i] = cnt;
8
        }
9
    }
10
    // 注意要把最后一个质数包括进去
11
    if (x > 1) pcnt[x] = 1;
12
    return pcnt;
13
}

注意我们可以直接在$n$上面做除法,而不用真的遍历到$\sqrt{n}$,因为发现一个质因数就可以把数据规模缩小,可以用哈希表把$log$消掉。这样,在最坏情况下这个算法时间复杂度是$O(\sqrt{n})$。当然,如果先把所有质数预处理出来,可以达到$O(\sqrt{\frac{n}{\log{n}}})$。

Pollard’s rho 算法

生日悖论 (Birthday Paradox)

在介绍Pollard’s rho算法之前,我想先介绍一个概率学中的小问题:一个房间中要有多少人,出现两个人生日相同的概率大于等于$50\%$?

如果你是第一次看到这个问题,那么答案可能会让你惊讶:只需要$23$个人。这不是玄学,而是通过严谨的数学推理得出的,实验也证明了这一点。那么这个结果是如何得到的呢?

假设一年有$365$天,且每个人的生日是独立且均匀分布的,那么第一个人在第$x$天生日的概率是$\frac{1}{365}$,我们记为$Pr\{b_1=x\}$。那么现在来了第二个人,它的生日在第$x$天的概率$Pr\{b_2=x\}$也是$\frac{1}{365}$,此时两个人都在第$x$天的概率是
$$
Pr\{b_1=x \textbf{ and } b_2=x\} = \frac{1}{365^2}
$$
但是这个$x$我们可以取$365$个,于是两个人生日都在同一天的概率就是
$$
Pr\{b_1=b_2\} = 365\times \frac{1}{365^2} = \frac{1}{365}
$$
知道了任意两人生日相同的概率,就可以扩展到$k$个人中的两人生日相同的期望值,由于$k$个人中选$2$个,所以有:
$$
E = {k\choose2} \frac{1}{365} = \frac{k(k-1)}{730}
$$
所以期望存在至少一对生日相同的两人只需要$k=28$,$50\%$概率的计算由于比较复杂就不放上来了,但是可以推测,这个值一定会比$28$少。

于是我们可以知道,生日碰撞的概率是$O(\sqrt{n})$级别的,而不是我们下意识的$O(n)$。

随机化算法

Pollard’s rho 算法的基本原理就是从待分解数$x$中任意找出两个数$a, b$,计算$gcd(|a - b|, x)$,当他们的$gcd$不是$1$的时候,我们就找到了$x$的一个因数。为什么是取两个呢?这就要利用我们上面说到的生日悖论了,因为只取一个$a$,此时$a|x$的几率很小,但是如果我们取两个数做差,出现$|a-b|$与$x$有公约数的几率就大很多了。事实上Pollard’s rho算法输出$x$的一个约数$p$的期望运行时间是$O(\sqrt{p})$。

那么如何选取这两个数呢?我们可以随机选取,但是Pollard’s rho算法使用的是一种二次剩余系的伪随机数生成函数
$$
f(x) = x^2 + c \pmod{n}
$$
为什么用这个函数呢?我没有在算法导论上看到具体说明,但是这个函数应该是对于Pollard’s rho算法的时间效率最优的函数。由于这个函数最后一定会出现循环,长得像希腊字母$\rho$,因此得名。


正因为这个性质,在出现环的时候我们也需要及时的退出并且换下一个随机数。我们可以使用Floyd判圈算法来判断是否出现了环,大意就是两个指针一个走一步一个走两步,如果有环那么两个指针一定会碰上,于是我们就可以有代码:

1
inline ll f(ll x, ll c, ll p) { return (modmul(x, x, p) + c) % p; }
2
ll pollard_rho(ll x) {
3
    ll c = mt() % (x - 1) + 1;
4
    ll l = f(mt() % x, c, x), r = f(l, c, x);
5
    // l是慢的指针,r是快的
6
    while (l != r) {
7
        ll g = __gcd(abs(l - r), x);
8
        if (g > 1) return g;
9
        l = f(l, c, x), r = f(f(r, c, x), c, x);
10
    }
11
    return x;
12
}

如果$x$是个质数,那么我们可以用上一章的Miller–Rabin算法提前判掉,这样整个算法的期望时间复杂度就是$O(n^{0.25}\log{n})$了。

常数优化

上面的写法虽然时间复杂度是正确的,但是对于这道题还是不够的:
例题:【模板】Pollard-Rho算法

我们还需要一些常数优化,Pollard’s rho算法的瓶颈就在于需要多少次才能得到一个约数$p$,以及找不到约数的时候做$gcd$的时间。前者并不好优化(玄学),但是减少$gcd$的次数却是可以做到的。

一个普遍的做法就是利用倍增的思想,先走完$2^k$步以后把$|a-b|$乘积统一进行$gcd$,然后换一个起点继续增加步数。如果出现乘积为$0$就说明我们遇到环了,此时退出并寻找下一个随机数。

1
ll pollard_rho(ll x) {
2
    ll c = (ll)mt() % (x - 1) + 1;
3
    ll s = 0, t = (ll)mt() % (x), val = 1;
4
    // i是步长
5
    for (int i = 1;; i <<= 1, s = t, val = 1) {
6
        for (int z = 1; z <= i; z++) {
7
            t = f(t, c, x);
8
            val = modmul(val, abs(t - s), x);
9
            // 如果遇到环了
10
            if (!val) return x;
11
            // 为了避免步数过长导致出不来结果,所以手动设定每127步计算一次
12
            if (!(z % 127)) {
13
                ll g = __gcd(val, x);
14
                if (g > 1) return g;
15
            }
16
        }
17
        ll g = __gcd(val, x);
18
        if (g > 1) return g;
19
    }
20
    return x;
21
}