模n的平方根算法

根存在性与数量

计算模 $n$ 的平方根其实就是求解方程
$$
x^2=y\ mod\ n\ ,\ gcd(y,n)=1
$$
那么显然，方程要有解，$y$ 必须是模 $p$ 的二次剩余。

当 $n$ 是素数时，根据欧拉判别法，$y$ 需要满足：
$$
y^{\frac{p-1}{2}}=1\ mod\ p
$$
否则无解。

对于 $n$ 是素数的情况，方程有 0 或 2 个解

对于 $n$ 是奇合数的情况，方程有 0 或 $2^l$ 个解（ $n=\sum_{i=1}^lp_i^{e_i}$ ）

$p=3\ mod\ 4$

关于这类素数，有一个通用公式
$$
x=\pm y^{\frac{p+1}{4}}\ mod\ p
$$
这个很好理解。~~其实我不会证明，但是可以用反证法（指反过来证明），计算一下 $x^2$ 发现正好是 $y$ 。~~

$p=1\ mod\ p$

这类素数的公式其实也是上面那个，但是问题是当 $p=1\ mod\ p$ 时，$\frac{p+1}{4}$ 并不是整数，而我们无法在模 $p$ 的情况下直接计算分数次幂（不然这篇东西就不必要存在了）。

于是我们有了 $Cipolla$ 算法。~~（下面是魔法）~~

首先我们对方程进行亿些变形
$$
x=\sqrt y=(a^2-(a^2-y))^{1/2}\ mod\ p
$$
令 $w^2 =a^2-y\ mod\ p$ ，并且 $(a^2-y)^{\frac{p-1}{2}}=-1\ mod\ p$ ，得到
$$
x=(a^2-w^2)^{1/2}=((a+w)(a-w))^{1/2}\ mod\ p
$$
又有
$$
w^p=w^{p-1}*w=(a^2-n)^{\frac{p-1}{2}}*w=-w\ mod\ p
$$
$$
a^p=a^{p-1}*a=a\ mod\ p
$$

所以
$$
x=((a^p+w^p)(a+w))^{1/2}\ mod\ p
$$
还有 $(a+w)^p=a^p+w^p\ mod\ p$ ，故
$$
x=(a+w)^{\frac{p+1}{2}}\ mod\ p
$$
现在我们的指数已经变成整数了，但问题是能不能求出这样的 $a$ 、$w$ ？

先来看 $a$ ， $a$ 是由 $(a^2-y)^{\frac{p-1}{2}}=-1\ mod\ p$ 决定的，换句话说， $a^2-y$ 是 $p$ 的二次非剩余。这样的 $a$ 是很好找的，因为模 $p$ 的既约剩余系中有一半是它的二次非剩余，我们只要随机找一个 $a$ ，然后带入检验即可，这一步不会花费太久。

然后来看 $w$ ，根据前面提到的， $w$ 是二次非剩余 $a^2-y$ 的平方根，显然 $w$ 不存在。另一方面， $w^2$ 是确确实实存在的。那我们就算不下去了吗？显然不是，我们再回头看看 $x=(a+w)^{\frac{p+1}{2}}\ mod\ p$ ，如果我们对其进行展开，会得到 $c+wd\ ,\ c、d\in R$ 这样的一个数（是不是很像复数）。按道理 $w$ 是不存在的，所以 $c+wd$ 也不应该存在，但事实是 $x$ 是存在的（否则相当于方程无解），那么答案只有一个——虚部为零，即 $d=0$ 。

所以我们得出结论，计算结果与$w$ 无关。

至于结果的计算方法，我们可以对它直接按照二项式定理计算。当然了，还有一种比较有趣的算法。前面说了 $c+dw$ 的形式很像复数，那我们也可以按照复数的计算来做。只不过这样我们需要自定义一个“复数”乘法的规则（毕竟计算机不知道），由于还涉及到了乘方，所以我们还需要定义一下快速幂算法。综合来看两种计算过程应该差不多，后者可以将更多复杂的计算过程交给计算机处理，个人感觉好一点。

下面放上脚本

#x^2=c mod p
def Cipolla_algorithm(c,p):
    # 定义模p复数域上的乘法
    def multiplication(x1,y1,x2,y2,w,p):
        x=(x1*x2+y1*y2*w)%p
        y=(x1*y2+x2*y1)%p
        return x,y

    # 获取w,使得w=-1mod(p),w是复数元的平方
    def get_w(c,p):
        a=randint(1,p)
        w=pow(a, 2) - c
        while pow(w,(p - 1)//2,p)!=p-1:
            a=randint(1,p)
            w=pow(a,2)-c
        return w,a
    #主体部分
    w,a=get_w(c,p)
    power=(p+1)//2
    x1=a
    y1=1
    x=1
    y=0
    while(power!=0):
        if(power!=(p+1)//2):
            x1,y1=multiplication(x1,y1,x1,y1,w,p)
        if(power & 1):
            x,y=multiplication(x,y,x1,y1,w,p)
        power>>=1
    return x

$n=pq$

这里拿它做一个示例，对于多因子 $n$ 的解法也差不多的。

首先根据中国剩余定理，方程 $x^2=y\ mod\ n\ ,\ n=pq$ 可以写成下面方程组
$$
\begin{cases}x^2=y\ mod\ p\\
x^2=y\ mod\ q\end{cases}
$$

然后分别解出来，我们得到了2组（4个）解(对于有 $l$ 种因子的 $n$ ，得到 $l$ 组)，按照中国剩余定理的操作，我们需要从每一组解中选出一个来进行最后的运算，所有最后一共可以得到 $2^l$ 个解。

不难看出，只有中国剩余定理的每个方程都有解，原方程才有解。