生日悖论

假设一年有$n$天，房间中有$k$人，每个人的生日在这$n$天中，服从均匀分布，两个人的生日相互独立
问至少要有多少人，才能使其中两个人生日相同的概率达到$p$

解：考虑$k\le n$
设$k$个人生日互不相同为事件$A$，则
$$P\left(A\right)=\frac{n}{n}\frac{n-1}{n}\cdots \frac{n-k+1}{n}$$
由题意$P\left(A\right)\le 1-p$,由$1+x\le e^x$
$$P\left(A\right)\le e^{-\frac{1}{n}}{e}^{-\frac{2}{n}}\cdots e^{\frac{k-1}{n}}=e^{-\frac{k\left(k-1\right)}{2n}}\le 1-p$$
解得
$$k\ge \frac{1+\sqrt{1-2n\ln \left(1-p\right)}}{2}$$

当$n=365,p=\frac{1}{2}$时，$k=23$,也就是说一个房间中至少23人就能使其中两个人生日相同的概率达到$50\%$
由于这个事实十分反直觉，故称为一个悖论

Pollard-Rho算法

Pollard-Rho算法是一种用于快速分解非平凡因数的算法
通过$f\left(x\right)=\left(x^2+c\right)\mathrm{mod}n$来生成一个随机数序列 { x i } \left\{x_i\right\} {xi​},其中$c$是一个随机数
随机取一个$x_1$,令$x_2=f\left(x_1\right),\cdots ,x_i=f\left(x_{i-1}\right)$,
这样产生的序列会形成一个$\rho$,也就是说会产生一个环

$\forall k\in {\mathbb{N}}_{+},gcd\left(k,n\right)\mid n$,只要选取适当的$k$使得$1<gcd\left(k,n\right)<n$,就能得到一个约数$gcd\left(k,n\right)$
满足这样的条件的$k$不少，$k$有若干质因子，每个质因子及其倍数都是可行的

由生日悖论，伪随机数序列中不同值的数量约为$O\left(\sqrt{n}\right)$(怎么算出来的其实我也不知道 )
设$m$为$n$的最小非平凡因子，显然$m\le \sqrt{n}$,令$y_i=x_i\mathrm{mod}m$
若$1\le c<n$,则
$$\begin{array}{rl}{y}_{i+1}& =x_{i+1}\mathrm{mod}m\\ & =\left(\left(x_i^2+c\right)\mathrm{mod}n\right)\mathrm{mod}m\\ & =\left(x_i^2+c\right)\mathrm{mod}m\\ & =\left({\left(x_i\mathrm{mod}m\right)}^2+c\right)\mathrm{mod}m\\ & =\left(y_i^2+c\right)\mathrm{mod}m\end{array}$$
于是我们可以得到新的序列 { y i } \left\{y_i\right\} {yi​}并且根据生日悖论可以得知序列中不同值的个数约为$O\left(\sqrt{m}\right)\le O\left(n^{\frac{1}{4}}\right)$

假设存在两个位置$i,j$，使得$x_i\ne x_j$且$y_i=y_j$,则$n\nmid \left|x_i-x_i\right|$且$m|\left|x_i-x_j\right|$
进而$1<gcd\left(\left|x_i-x_j\right|,n\right)<n$
因此我们可以通过$gcd\left(\left|x_i-x_j\right|,n\right)$获得$n$的一个非平凡因子

floyd判环

$a=f(a)$
$b=f(f(b)$
如果$a=b$则有环，就直接返回，更换个$c$再来一遍
如果$gcd\left(\left|a-b\right|\right)>1$,则得到了一个因数

洛谷P4718
__int128 + 优化gcd + 7个数字判断质数 + O2才能过（不开O2会T第13个点）

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <ctime>
using namespace std;
typedef long long LL;
LL gcd(LL x, LL y){
    if (!x) return y;
    if (!y) return x;
    LL t = __builtin_ctzll(x | y);
    x >>= __builtin_ctzll(x);
    do
    {
        y >>= __builtin_ctzll(y);
        if (x > y) swap(x, y);
        y -= x;
    } while (y);
    return x << t;
}
LL quick_pow(LL a, LL b, LL p) {
    LL res = 1;
    while (b) {
        if (b & 1)res = (__int128)res * a % p;
        a = (__int128)a * a % p;
        b >>= 1;
    }
    return res;
}
const LL bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
bool Rabin_Miller(LL n) {
    if (n < 3 || (n & 1) == 0)return n == 2;

    LL u = n - 1, t = 0;
    // n - 1 = u * (2^t)
    while ((u & 1) == 0) {
        u >>= 1;
        ++t;
    }
    for (int i = 0; i < 7; ++i) {
        LL a = bases[i] % n;
        if (a == 0)continue;
        LL v = quick_pow(a, u, n);
        if (v == 1)continue;
        LL s;
        for (s = 0; s < t; ++s) {
            //a^{u * 2^s}= -1 (mod n)
            if (v == n - 1)break;
            v = (__int128)v * v % n;
        }
        if (s == t)return false;
    }
    return true;
}
LL max_factor;
LL f(LL x, LL c, LL p) {
    return ((__int128)x * x % p + c) % p;
}
LL getRandom(const LL& a, const LL& b) {
    return ((1LL * rand() << 32LL) + 1LL * rand()) % (b - a + 1LL) + a;
}
LL Pollard_Rho_floyd(LL x) {
    if (x == 4)return 2;
    LL c = getRandom(3, x - 1);//[3,x-1]
    LL s = getRandom(0, x - 1);//[0,x-1]
    s = f(s, c, x);
    LL t = f(s, c, x);
    while (s != t) {
        LL d = gcd(abs(t - s), x);
        if (d > 1)return d;
        s = f(s, c, x);
        t = f(f(t, c, x), c, x);
    }
    return x;
}
void fac(LL x) {
    if (x <= max_factor || x < 2)return;
    if (Rabin_Miller(x)) {// x是质数
        if (x > max_factor) {
            max_factor = x;
        }
        return;
    }
    LL p = x;
    while (p >= x)p = Pollard_Rho_floyd(x);//找一个因子p
    while (x % p == 0)x /= p;
    fac(x);
    fac(p);
}
int main() {
    srand((unsigned)time(NULL));
    int T;
    scanf("%d", &T);
    while (T--) {
        max_factor = 1;
        LL n;
        scanf("%lld", &n);
        fac(n);
        if (n == max_factor)printf("Prime\n");
        else printf("%lld\n", max_factor);
    }
    return 0;
}

倍增优化

由于频繁使用gcd会导致复杂度上去，
我们考虑累积几次再算gcd
下面是$min\left(2^k-1,128\right)$次算一个gcd

洛谷P4718
__int128

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <ctime>
using namespace std;
typedef long long LL;
LL gcd(LL a, LL b){
    LL c;
    while (b){
        c = a % b;
        a = b;
        b =c;
    }
    return a;
}
LL quick_pow(LL a, LL b, LL p) {
    LL res = 1;
    while (b) {
        if (b & 1)res = (__int128)res * a % p;
        a = (__int128)a * a % p;
        b >>= 1;
    }
    return res;
}
const LL bases[] = {2, 325, 9375, 28178, 450775, 9780504, 1795265022};
bool Rabin_Miller(LL n) {
    if (n < 3 || (n & 1) == 0)return n == 2;

    LL u = n - 1, t = 0;
    // n - 1 = u * (2^t)
    while ((u & 1) == 0) {
        u >>= 1;
        ++t;
    }
    for (int i = 0; i < 7; ++i) {
        LL a = bases[i] % n;
        if (a == 0)continue;
        LL v = quick_pow(a, u, n);
        if (v == 1)continue;
        LL s;
        for (s = 0; s < t; ++s) {
            //a^{u * 2^s}= -1 (mod n)
            if (v == n - 1)break;
            v = (__int128)v * v % n;
        }
        if (s == t)return false;
    }
    return true;
}
LL max_factor;
LL f(LL x, LL c, LL p) {
    return ((__int128)x * x % p + c) % p;
}
LL getRandom(const LL& a, const LL& b) {
    return ((1LL * rand() << 32LL) + 1LL * rand()) % (b - a + 1LL) + a;
}
LL Pollard_Rho(LL x) {
    if (x == 4)return 2;
    LL c = getRandom(3, x - 1);//[3,x-1]
    LL s = getRandom(0, x - 1);//[0,x-1]
    s = f(s, c, x);
    LL t = f(s, c, x);
    for (int lim = 1; s != t; lim = std::min(128, lim << 1)) {
        LL val = 1;
        for (int i = 0; i < lim; ++i) {
            LL temp = (__int128)val * abs(s-t) % x;
            if (temp == 0)break;
            val = temp;
            s = f(s, c, x);
            t = f(f(t, c, x), c, x);
        }
        LL d = gcd(val, x);
        if (d > 1)return d;
    }
    return x;
}
void fac(LL x) {
    if (x <= max_factor || x < 2)return;
    if (Rabin_Miller(x)) {// x是质数
        if (x > max_factor) {
            max_factor = x;
        }
        return;
    }
    LL p = x;
    while (p >= x)p = Pollard_Rho(x);//找一个因子p
    while (x % p == 0)x /= p;
    fac(x);
    fac(p);
}
int main() {
    srand((unsigned)time(NULL));
    int T;
    scanf("%d", &T);
    while (T--) {
        max_factor = 1;
        LL n;
        scanf("%lld", &n);
        fac(n);
        if (n == max_factor)printf("Prime\n");
        else printf("%lld\n", max_factor);
    }
    return 0;
}