问题

求$x^b\equiv a\pmod n$的所有根，其中$\gcd(a,n)=1$且根的数目合适，$n$的分解$\prod\nolimits_i p_i^{e_i}$已知/能在合适时间内给出。

思路

对所有$p_i\ne 2$，令：

$r=\gcd(b,p_i-1)$

那么：

$\begin{align*}x^r&\equiv a^{(b/r)^{-1}\bmod (p_i-1)/r}\\&\equiv \delta\pmod {p_i}\end{align*}$

若$r\ne 1$，接着用AMM算法求$x^r\equiv \delta\pmod {p_i}$的根。

最后lift一下；

对于$p_i=2$，我将在后面单独讨论。

最最后用CRT合并。

AMM

接下来讨论$p\ne 2$的情形。AMM算法可用于求有限域中的平方根（即$r=2$），这里给出了对$r$的拓展。

论文只讨论了$p-1=r^ts,\gcd(r,s)=1$的情形，而实际上只要分解$r=\prod\nolimits_i p_i^{e_i}$，再对所有$p_i$使用$e_i$次该算法，就可以让算法适用于所有大小合适的$r$（不过这是效率最低的办法，每次都要重新算幂）。

那么下面假定$r$为素数，$p-1=r^ts,\gcd(r,s)=1$，我们有：

$\begin{align*}(\delta^{r^{-1}\bmod s})^r&\equiv\delta^{1+ks}\\&\equiv\delta^{1+k(p-1)/r^{t}}\pmod p\end{align*}$

只要求出：

$c\equiv(\delta^{1-r(r^{-1}\bmod s)})^{1/r}\pmod p$

就有：

$\begin{align*}x&\equiv\delta^{1/r}\\&\equiv \omega_r c\delta^{r^{-1}\bmod s}\pmod p\end{align*}$

其中$\omega_r$为任意一个$r$次单位根，而$\delta^{1-r(r^{-1}\bmod s)}$为某个$r^{t-1}$次单位根。

接着读论文里的算法，我突然发现如下事实：

任选$r$次非剩余$\rho$，$g\equiv \rho^{(p-1)/{r^{t-1}}}\pmod p$在乘法作用下生成一个阶为$r^{t-1}$的子群$G$，$G$的元素就是所有$r^{t-1}$次单位根。搁$G$上求离散对数不就完事儿了吗：

令：

$h=\log_g(\delta^{1-r(r^{-1}\bmod s)})$

则：

$c\equiv \rho^{h(p-1)/r^t}\pmod p$

随机选取$\rho$，$\rho$为$r$次非剩余，即$\rho^{(p-1)/r}\not\equiv 1\pmod p$的概率为$1-\frac{1}{r}$，这意味着期望在常数次内能选出一个$r$次非剩余$\rho$。

然后调各种解DLP的算法就好了。

为啥我要$r$是素数？论文中没有要求$r$为素数，从而产生了一个bug：

$r$非素时，虽然$\rho$为非$r$次剩余，但$\rho^{(p-1)/r}$生成的群的阶并不一定为$r$，这导致算法Step 4可能无法继续。一个办法是分解$r$，对所有素因子$p_i$检查$\rho^{(p-1)/{p_i}}$。~~不过我这里就直接要求$r$是素数了。~~

如果直接调Pohlig-Hellman+BSGS和最慢的整数乘法，且$r=\prod\nolimits_i p_i^{e_i}$，那么时间复杂度为（没包括选$\rho$）：

$\mathcal O(t\log^2 p\sum\nolimits_i e_i(t\log r+\sqrt {p_i}))$

Hensel Lifting

设：

$x_k=\delta^{1/r}\bmod {p^{2^k}}$

令：

$x_{k+1}=x_k+p^{2^k}y_k,\quad 0\le y_k<p^{2^k}$

那么：

$\begin{align*}x_{k+1}^r&\equiv (x_k+p^{2^k}y_k)^r\\&\equiv x_k^r+rx_k^{r-1}p^{2^k}y_k\\&\equiv \delta\pmod {p^{2^{k+1}}}\end{align*}$

从而：

$p^{2^k}y_k= ((rx_k^{r-1})^{-1}(\delta-x_k^r))\bmod p^{2^{k+1}}$

实际上不需要计算$(rx_k^{r-1})^{-1}\bmod {p^{2^{k+1}}}$，因为：

$\begin{align*}y_k&= (rx_k^{r-1})^{-1}\frac{(\delta-x_k^r)}{p^{2^k}}\bmod {p^{2^k}}\\p^{2^k}y_k&= (((x_k(r\delta)^{-1})\bmod {p^{2^k}})(\delta -x_k^r))\bmod p^{2^{k+1}}\end{align*}$

算出$x_{\lceil\log_2 e\rceil}\bmod p^e$就完事了。

然后算逆元也能说道说道，设：

$x_k=m^{-1}\bmod p^{2^k}$

如法炮制，可以得到：

$x_{k+1}=(x_k(2-mx_k))\bmod p^{2^{k+1}}$

用lift的话，总的来说就只是给复杂度后边加上个（假定用最慢的整数乘法）$\mathcal O(e^2\log r\log^2 p+\log^3 p)$，不然得给前面那部分乘个$e^2$，就有点难受了。

“Even” Hard

关于$p=2$该咋整，我除了在Hensel基础上加点搜索以外真想不到别的法子。

不过Sage这里有一个没搜索的实现，注释说是”a variant of Hensel lifting”，不过我暂时没找到这个”variant”出处在哪，这源代码那个for也没能看懂。

等我懂了再补完吧。