Loading [MathJax]/jax/element/mml/optable/BasicLatin.js

二次剩余Cipolla算法学习笔记

 


对于同余式

x2n(modp)

若对于给定的n,P,存在x满足上面的式子,则乘n在模p意义下是二次剩余,否则为非二次剩余

我们需要计算的是在给定范围内所有满足条件的x,同时为了方便,我们只讨论p是奇质数的情况

前置定理

  • x2(x+p)2(modp)

证明:x2x2+2xp+p2(modp)显然成立

  • 对于x2n(modp),除n=0外,总共有p12n使得该方程有解

我局的参考资料里对于这条性质的证明漏洞很大,所以下面的是自己yy的

根据第一个前置定理的式子,我们只需讨论x[1,p1]即可(当x=0时对应了n=0的特殊情况)

一个显然的性质是

x2(px)2(modp)

那么当x[1,p12]我们可以取到所有解。

接下来我们只需要证明当x[1,p12]x2mod均两两不同

可以用反证法,若存在不同的u, v满足u^2 \equiv v^2 \pmod p

那么有(u + v)(u - v) \equiv 0 \pmod p

显然-p < u + v < p-p < u - v < pu + v \not = 0, u - v \not = 0,故该假设不成立,故原命题成立。

Q.E.D

  • 勒让德符号(Legender symbol)

(\frac{a}{p}) = \begin{cases} 1 , &\text{a在模$p$意义下是二次剩余}\\ -1, &\text{a在模$p$意义下是非二次剩余}\\ 0, &\text{a mod p = 0} \end{cases}

这个东西的分布大概是这个样子

image

计算公式

我局的这个公式就是构造出来的

(\frac{a}{p}) = a^{\frac{p - 1}{2}} \pmod p

证明:

费马小定理:对于任意互质的x, p,有x^{p - 1} = 1 \pmod p

一条同余式的性质:若a^k \equiv b^k \pmod p,那么a^{kx} \equiv b^{xk} \pmod p

然后直接把这玩意儿带到x^2 \equiv a \pmod p里就行了

这里简单的写一下:

首先要明确我们的目的,我们现在要验证这个公式的正确性,也就是说我们要证明当a^{\frac{p-1}{2}}=1 \pmod p时满足条件的x存在,当a^{\frac{p-1}{2}}= -1 \pmod px不存在,当a^{\frac{p-1}{2}}= 0 \pmod pa\mod p = 0

  1. a^{\frac{p-1}{2}}=1 \pmod p

我们假设有x^2 \equiv a \pmod p

x^{2\frac{p-1}{2}} \equiv a^{\frac{p-1}{2}} \pmod p

x^{p-1} \equiv 1 \pmod p

根据费马小定理x显然存在,因此a是模p意义下的二次剩余

  1. a^{\frac{p-1}{2}}= -1 \pmod p

假设有x^2 \equiv a \pmod p

同理可知

x^{p-1} \equiv -1 \pmod p

显然x不存在,因此a不是模p意义下的二次剩余

  1. a^{\frac{p-1}{2}}= 0 \pmod p

显然有a \bmod p = 0

Cipolla算法

算法流程

这个算法其实用两句话就能说完,但是背后的理论却非常高深(对于我这种菜鸡而言)。

  1. 首先使用随机的方法找到一个(\frac{a^2 - n}{p}) = -1,记\omega = \sqrt{a^2-n}

  2. 那么x \equiv (a + w)^{\frac{p+1}{2}} \pmod p

做完了。。。期望复杂度O(\log^2 n)

但是实际上实现起来并没有这么简单,因为要自定义类似于虚数的乘法/幂运算

算法理论

首先要有一点抽代基础(群/环/域什么的要知道定义)

我们来逐步分析这个算法(按照我的叙述风格应该是从发明者的角度出发一步一步推出这玩意儿来,但是十分抱歉我实在是搞不明白他当时的脑回路qwq)

对于第一步,根据前面的定理,如果在[1, p]内随机,每次有\frac{1}{p}*\frac{p-1}{2}的概率找到一个解,那么期望步数大约为两次,因此复杂度是可以保证的。

但是找到这个东西有什么用呢?。如果我们把之前的数域记做\mathbf F_p\omega在这个数域下是不能开根的,但是我们可以构造一个新的数域\mathbf F_p,使得\omega\mathbf F_{p2}下能够开根。类比于-1在复数域下能够表示为\sqrt{-1}一样。

这样的话\mathbf F_{p2}内的数都可以写作a + k\omega的形式。可以证明这玩意儿确实是个合法的域,证明过程,同时也可以证明在\mathbf F_{p2}下得到的解在\mathbf F_{p1}下也成立,同时最后的答案中\omega的系数一定为0

现在来简单说明一下为什么x \equiv (a+\omega)^{\frac{p+1}{2}}

先来了解两个性质

  • \omega^p \equiv -\omega \pmod p

证明:

\begin{aligned} \omega^p &= (a^2-n)^{\frac{p}{2}}\\ &= (a^2 - n)^{\frac{p - 1}{2}} (a^2 - n)^{\frac{1}{2}}\\ &= -\omega \end{aligned}

  • (a + b)^p \equiv a^p + b^p \pmod p

证明就直接考虑二项式定理中的组合数展开,发现除了第一项和最后一项之外都无法把n!消掉。

那么要证明x \equiv (a+\omega)^{\frac{p+1}{2}},实际上我们只需要证明(a+\omega)^{p+1}\equiv n \pmod p就行了

\begin{aligned} &(a + \omega)^{p + 1}\\ =&(a + \omega)^p(a + \omega)\\ =&(a - \omega)(a + \omega)(\text{根据费马小定理$a^p \equiv p \pmod p$})\\ =&(a^2 - \omega^2)\\ =&(a^2 - (a^2 - n))\\ =&n \end{aligned}

算法的大概思想就讲完了,下面煮个栗子~。

对于x^2 \equiv n \pmod p

假设此时p=13, n = 10

首先要找到一个a满足(\frac{a^ - 10}{13}) = -1,然后脸黑的attack在经过1e9 +7次尝试后终于找到了一个a =2它满足条件,因为(\frac{7}{13}) = -1此时\omega = \sqrt{a^2 - n} = \sqrt{-6}

按照老祖宗讲给我们的

x \equiv (2 + \sqrt{-6})^{7} \pmod {13}

\begin{aligned} &\left(2+{\sqrt {-6}}\right)^{2}=4+4{\sqrt {-6}}-6=-2+4{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{4}=\left(-2+4{\sqrt {-6}}\right)^{2}=-1-3{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{6}=\left(-2+4{\sqrt {-6}}\right)\left(-1-3{\sqrt {-6}}\right)=9+2{\sqrt {-6}}\\ &\left(2+{\sqrt {-6}}\right)^{7}=\left(9+2{\sqrt {-6}}\right)\left(2+{\sqrt {-6}}\right)=6. \end{aligned}

然后不难发现36 \equiv 10 \pmod {13}

同时因为平方的性质,-x也是一个合法解,因此-6 + 13 = 7也是合法的

最后有一个小问题就是为什么最后\omega的系数一定是0,参考资料中给出的解释我实在是不能理解,如果有看得懂的大佬欢迎给本菜鸡讲一下qwq

image

代码模板

#include<bits/stdc++.h>  
using namespace std;
const int mod = 13;
namespace TwoRemain {
template <typename A, typename B> inline int add(A x, B y) {if(x + y < 0) return x + y + mod; return x + y >= mod ? x + y - mod : x + y;}
template <typename A, typename B> inline void add2(A &x, B y) {if(x + y < 0) x = x + y + mod; else x = (x + y >= mod ? x + y - mod : x + y);}
template <typename A, typename B> inline int mul(A x, B y) {return 1ll * x * y % mod;}
template <typename A, typename B> inline void mul2(A &x, B y) {x = (1ll * x * y % mod + mod) % mod;}
int fmul(int a, int p, int Mod = mod) {
    int base = 0;
    while(p) {
        if(p & 1) base = (base + a) % Mod;
        a = (a + a) % Mod; p >>= 1;
    }
    return base;
}
int fp(int a, int p, int Mod = mod) {
    int base = 1;
    while(p) {
        if(p & 1) base = fmul(base, a, Mod);
        p >>= 1; a = fmul(a, a, Mod);
    }
    return base;
}
int f(int x) {
    return fp(x, (mod - 1) >> 1);
}
struct MyComplex {
    int a, b;
     int cn;
    MyComplex operator * (const MyComplex &rhs)  {
        return {
            add(fmul(a, rhs.a), fmul(cn, fmul(b, rhs.b, mod))),
            add(fmul(a, rhs.b), fmul(b, rhs.a)),
            cn
        };
    }
};
MyComplex fp(MyComplex a, int p) {
    MyComplex base = {1, 0, a.cn};
    while(p) {
        if(p & 1) base = base * a;
        a = a * a; p >>= 1;
    }
    return base;
}
int TwoSqrt(int n) {
    if(f(n) == mod - 1) return -1;
    if(f(n) ==  0) return  0;
    int a = -1, val = -1;
    while(val == -1) {
        a = rand() << 15 | rand();
        val = add(mul(a, a), -n);
        if(f(val) != mod - 1) val = -1;
    }
    return fp({a, 1, val}, (mod + 1) / 2).a;
}
}
using namespace TwoRemain;
signed main() {
    cout << TwoSqrt(10);
    return 0;
}

参考资料

二次剩余Cipolla算法学习小记

Legendre symbol

posted @   自为风月马前卒  阅读(1843)  评论(8编辑  收藏  举报
编辑推荐:
· 深入理解 Mybatis 分库分表执行原理
· 如何打造一个高并发系统?
· .NET Core GC压缩(compact_phase)底层原理浅谈
· 现代计算机视觉入门之:什么是图片特征编码
· .NET 9 new features-C#13新的锁类型和语义
阅读排行:
· Spring AI + Ollama 实现 deepseek-r1 的API服务和调用
· 《HelloGitHub》第 106 期
· 数据库服务器 SQL Server 版本升级公告
· 深入理解Mybatis分库分表执行原理
· 使用 Dify + LLM 构建精确任务处理应用
历史上的今天:
2018-03-27 利用MingW检验程序运行内存
2017-03-27 2821 天使之城
2017-03-27 3139 栈练习3
2017-03-27 3138 栈练习2
2017-03-27 3137 栈练习1
2017-03-27 中缀表达式值
2017-03-27 车厢调度(train.cpp)

Contact with me

点击右上角即可分享
微信分享提示