这篇题解不使用 Pollard-rho 算法。
下面介绍二次筛法（Quadratic Sieve）。它的时间复杂度约为 O(e^{3c\sqrt{\ln n\ln\ln n}})，这里 c 是一个常数，理论上当 n 足够大时应该取 c=\dfrac{\sqrt{2}}{4}。在本题数据范围内取 c=0.46 最好。
参考资料
前置知识：
线性筛
Euler 判别法和 Cipolla 算法
Guass-Jordan 消元法
Miller-Rabin 素性测试
二次筛法主要用于分解 10^{30} 以上的大整数，在本题数据范围内比 Pollard-rho 慢一些。
1 算法思想
1.1 Fermat 分解法
给定正整数 n，设 x^2-y^2=n，则有 (x+y)\mid n。于是可以从 y=\lceil\sqrt{n}\rceil 开始逐个尝试 y，直到得到使得 x=\sqrt{n+y^2}\in\mathbb{N} 的 y。这个方法对于因子很接近 \sqrt{n} 的 n 有效。
1.2 二次筛法流程
尝试推广 Fermat 分解法：如果能找到 x^2\equiv y^2\pmod n，那么 \gcd(x+y,n) 大概率是 n 的一个因子。于是只需找到足够多的 x 和 y 满足 x^2\equiv y^2\pmod n 就很可能分解出 n。
考虑函数 Q(x)=(x+\lfloor\sqrt{n}\rfloor)^2-n。记 \tilde{x}=(x+\lfloor\sqrt{n}\rfloor)，有 Q(x)\equiv\tilde{x}^2\pmod n。如果已经找到 x_1,x_2,\dots,x_k 和 y 使 \left(\prod_{i=1}^kQ(x_i)\right)=y^2，那么 \left(\prod_{i=1}^k\tilde{x_i}\right)^2\equiv y^2\pmod n，于是 \gcd(\left(\prod_{i=1}^k\tilde{x_i}\right)+y,n) 就可能是 n 的一个因子。
为了找到 x_1,x_2,\dots,x_k 使 \prod_{i=1}^kQ(x_i) 为平方数，首先给定一组小素数（它们称为因子基）p_1,p_2,\dots,p_B，然后考虑那些在因子基上光滑（即素因子都在因子基内）的 Q(x)。设 Q(x_i)=\prod_{j=1}^Bp_j^{\alpha_{i,j}}，则只需满足对任意 j 有 \left(\sum_{i=1}^k\alpha_{i,j}\right)\equiv0\pmod2 就能保证 \prod_{i=1}^kQ(x_i) 是平方数。
由此得到二次筛法：
-
选定一组因子基。
-
得到足够多在因子基上光滑的 Q(x_i)，写出 \alpha_{i,j}。
-
选出一些 i_m，使得 \left(\sum_{m=1}^k\alpha_{i_m,j}\right)\equiv0\pmod2 对于任意 j 成立。
-
根据这些 i_m 计算出 \gcd(\left(\prod_{m=1}^k\tilde{x}_{i_m}\right)+\sqrt{\prod_{m=1}^kQ(x_{i_m})},n)。
2 算法细节
2.1 因子基的选取与一些预处理
这一步要选取合适的 p_1,p_2,\cdots,p_B，并准备进行二次筛需要的预处理。
因子基的大小 B\approx e^{c\sqrt{\ln n\ln\ln n}}，其中 c 的值见本文第二段。
要让 p\mid Q(x) ，即 \tilde{x}^2\equiv n\pmod p，故 n 为模 p 的二次剩余。取 p_1=2，p_2,\cdots,p_B 为前 B-1 个满足条件的奇素数即可。这用线性筛和 Euler 判别法容易实现。
选取好因子基后，为了方便之后快速判断 p\mid Q(x) 是否成立，还需要计算出模 p 意义下 n 的平方根。这可以用 Cipolla 算法。另外还需要提前计算好 \ln p 的值，以便判断 Q(x) 是否光滑。
2.2 二次筛
这一步要筛出光滑的 Q(x)。
从 x=1 开始依次枚举 Q(x) 。对于每个 x，利用提前处理好的数据判断 p_j\mid Q(x) 对于每个 j 分别是否成立，如果成立则从 \ln Q(x) 中减去 \ln p_j 。如果最终减去的结果很接近 0 ，那么 Q(x) 就是光滑的。如果结果大于 0 ，也有可能 Q(x) 中含有不止一个 p_j ，但为了省时，除了含有多个 p_1=2 的情况被考虑以外，直接忽略含有多个其它 p_j 的情况。如果 Q(x_i) 为光滑的，记录每个 p_j 是否为 Q(x_i) 的因数，是则 \alpha_{i,j}=1，否则 \alpha_{i,j}=0。这一步在筛出 B+trt 个 Q(x) 后可以停止，其中 trt 代表后面可以尝试的 \gcd(\cdot,n) 的个数。trt 如果过小可能会分解不出来，过大会需要的时间更多，取 trt=15 比较合适。
程序运行的时间主要在这一步。可以证明大约需要枚举 O(B^3) 次。
2.3 消元
这一步要对 \alpha_{i,j} 矩阵进行初等变换（\alpha_{i,j} 是 \bmod 2 意义下的，所以这里的初等变换指某一行 \operatorname{xor} 上另一行），从
\begin{bmatrix}\alpha_{1,1}&\alpha_{2,1}&\alpha_{3,1}&\cdots&\alpha_{B+trt,1}\\\alpha_{1,2}&\alpha_{2,2}&\alpha_{3,2}&\cdots&\alpha_{B+trt,2}\\\alpha_{1,3}&\alpha_{2,3}&\alpha_{3,3}&\cdots&\alpha_{B+trt,3}\\\cdots&\cdots&\cdots&\cdots&\cdots\\\alpha_{1,B}&\alpha_{2,B}&\alpha_{3,B}&\cdots&\alpha_{B+trt,B}\end{bmatrix}
变为
\begin{bmatrix} 1&0&0&\cdots&0&\alpha'_{1,1}&\alpha'_{2,1}&\cdots&\alpha'_{trt,1}\\0&1&0&\cdots&0&\alpha'_{1,2}&\alpha'_{2,2}&\cdots&\alpha'_{trt,2}\\0&0&1&\cdots&0&\alpha'_{1,3}&\alpha'_{2,3}&\cdots&\alpha'_{trt,3}\\\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\0&0&0&\cdots&1&\alpha'_{1,B}&\alpha'_{2,B}&\cdots&\alpha'_{trt,B}\end{bmatrix}
的形式，这里每一列代表一个 Q(x)，每一行代表因子基里的一个素数。
对于每个 1\le q\le trt，选 i_1,i_2,\cdots,i_{k-1} 为满足 \alpha'_{q,i}=1 的所有 i，i_k=q。那么在消元后的矩阵中，\left(\sum_{m=1}^k\alpha_{i_m,j}\right)\equiv0\pmod2 对于任意 j 成立，于是这个式子对原矩阵也成立。
这一步可以用在域 \mathbb{F}_2 上的 Gauss-Jordan 消元法，复杂度与 2.3 节同样是 O(B^3)。
2.4 分解
这一步要用上面的消元结果得出 n 的一个非平凡因子。
对于每个 1\le q\le trt：
设
\alpha_j=\dfrac{\sum_{m=1}^k\alpha_{i_m,j}}{2}
X=\prod_{j=1}^Bp_j^{\alpha_j}
Y=\prod_{m=1}^k\tilde{x}_{i_m}=\prod_{m=1}^k(x_{i_m}+\lfloor\sqrt{n}\rfloor)
由于
\begin{aligned}X^2&\equiv \prod_{j=1}^Bp_j^{\sum_{m=1}^k\alpha_{i_m,j}}\\&\equiv\prod_{j=1}^B\prod_{m=1}^kp_j^{\alpha_{i_m,j}}\\&\equiv\prod_{m=1}^k\prod_{j=1}^Bp_j^{\alpha_{i_m,j}}\\&\equiv\prod_{m=1}^kQ(x_{i_m})\\&\equiv\prod_{m=1}^k\tilde{x}_{i_m}^2\\&\equiv\left(\prod_{m=1}^k\tilde{x}_{i_m}\right)^2\\&\equiv Y^2\pmod n\end{aligned}
所以 c=\gcd(X+Y,n) 可能是 n 的一个因数。
每次尝试有至少 \frac{1}{2} 的概率成功，所以分解成功的概率至少有 1-2^{-trt}=1-2^{-15}>0.9999。
3 代码实现
完整代码
需要开 O2。代码中有一个小优化可以使运行时间减半：同时筛 Q(x) 和 -Q(-x)，在 \alpha 矩阵中增加一行记录是否为负数，消元后取偶数个负数即可。