数论内容简要整理
部分内容引自清华大学秦岳《初等数论》,dengyaotriangle的博客,_WZT_ 的博客,niiick的博客,超十一维的博客
算术基本定理
任何一个大于1的自然数$N$,如果$N$不为质数,那么$N$可以唯一分解成有限个质数的乘积
$$N=P_1^{a_1}*P_2^{a_2}*P_3^{a_3}......P_n^{a_n}$$
这里$P1<P2<P3......<Pn$均为质数,其中指数$a_i$是正整数。这样的分解称为$N$的标准分解式
素数无限定理
正整数集中包含无限个素数
证明:假设素数有限,设其为$p_1...p_n$,构造$s=1+\prod p_i$,若s是素数,矛盾;若s是合数,则$p_1...p_n$都不是s的约数,与算术基本定理矛盾
Eratosthnes素数筛法
设置一从2开始的数表,若该数没有被划掉,则其为一个素数,将其所有倍数划掉
复杂度$O(nlogn)$
优化版埃氏筛关键部分代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
p[1]=1; //bool p[i]表示i是否为素数 int sq=sqrt(n); for(int i=2;i<=sq;i++) if(!p[i]) for(int j=i*i;j<=n;j+=i) p[j]=1;
欧拉线性筛法
维护每个数的最小素因子$mindiv$,从2开始,若其$mindiv$未知,则为素数,将其小于本身$mindiv$的素数倍的数设置成素数。每个数只会被其最小素因子筛一次。复杂度$O(n)$
线性筛关键部分代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
p[1]=1; for(int i=2;i<=n;i++) { if(!p[i]) pr[++cnt]=i; //pr为素数表 for(int j=1;j<=cnt&&pr[j]*i<=n;j++) { p[i*pr[j]]=1; //以pr[j]作为最小质因子更新p[i*pr[j]] if(i%pr[j]==0) break; /*i的mindiv比pr[j]更小,最小质因子不再是pr[j] ,若继续更新则会出现冗余操作*/ } }
利用线性筛实现的高效分解质因数
对于数$x$每次取其最小质因子,更新$x$为$x/mindiv$,当$x=1$时分解结束,时间复杂度$O(logn)$
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
s[1]=1; for(int i=2;i<=n;i++) { if(!s[i]) p[++cnt]=i,g[i]=cnt; for(int j=1;j<=cnt&&i*p[j]<=n;j++) { s[i*p[j]]=1,g[i*p[j]]=j; if(i%p[j]==0) break; } }
朴素分解质因数
$O(\sqrt n)$的分解质因数,对于单个数的分解质因数和积性函数单值求解具有重要作用
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
for(int i=2;i*i<=n;i++) while(n%i==0) //这里的i一定是素数,不妨利用反证法证明 //若i为合数,设i有素因子j小于i,j在前面的过程中已被筛去,n%j!=0故n%i!=0,与条件矛盾,故i必为素数 n/=i;
更相减损术
$gcd(a,b)=gcd(a,a-b)$
欧几里得算法
$gcd(a,b)=gcd(b,a\%b)$
裴蜀定理(贝祖定理)
若$a,b$是整数,且$(a,b)=d$,那么对于任意的整数$x,y$,$ax+by$都一定是$d$的倍数,特别地,一定存在整数$x,y$,使$ax+by=d$成立。
推论:$a,b$互质的充要条件是存在整数$x,y$使$ax+by=1$
推论:设$a_1,a_2,a_3......a_n$为$n$个整数,$d$是它们的最大公约数,那么存在整数$x_1......x_n$使得$x_1*a_1+x_2*a_2+...x_n*a_n=d$
扩展欧几里得算法(exgcd)
$exgcd$用于解决二元一次不定方程的解的相关问题
对于不定方程$ax+by=c$
由裴蜀定理可知$gcd(a,b)|(ax+by)$,故若$c$不是$gcd(a,b)$的倍数则无整数解
引理:$ax+by=gcd(a,b)$的特解求法
显然存在$x1,y1$使得$ax+by=bx_1+(a\%b)y_1$成立
由模运算的性质可知$a\%b=a-b*\left\lfloor\dfrac{a}{b}\right\rfloor$
可得$ax+by=bx_1+(a-b*\left\lfloor\dfrac{a}{b}\right\rfloor)y_1$
$=ay_1+b*(x_1-\left\lfloor\dfrac{a}{b}\right\rfloor y_1)$
则有$x=y_1,y=x_1-\left\lfloor\dfrac{a}{b}\right\rfloor y_1$
由欧几里得算法得迭代过程的终止状态为$a=gcd(a,b),b=0$
此时有$gcd(a,b)*x=gcd(a,b)$,取$x=1$即可,此时任取$y$均可保证该式成立,为避免溢出常取$y=0$
该过程的代码实现
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
void exgcd(int a,int b,int &x,int &y) { if(!b) { x=1,y=0; return; } exgcd(b,a%b,y,x); y-=x*(a/b); }
设$x_0,y_0$为$ax+by=gcd(a,b)$的一组特解,即有$ax_0+by_0=gcd(a,b)$
则有原方程的一组特解
$$a\frac{x_0c}{gcd(a,b)}+b\frac{y_0c}{gcd(a,b)}=c$$
$$x_1=\frac{x_0c}{gcd(a,b)},y_1=\frac{y_0c}{gcd(a,b)}$$
构造原方程的通解形式
易得对于任意常数$d$,有$a(x_1+db)+b(y_1-da)=c$成立。为保证通解均为整数解,只需使$db,da$均为整数即可;而为保证得到通解,$d$须取满足上述条件的最小值,而易得该最小值为$\frac{1}{gcd(a,b)}$
可得通解:$x=x_1+\frac{kb}{gcd(a,b)},y=y_1+\frac{ka}{gcd(a,b)}$其中$k$为任意整数
并有性质:随$k$增大,$x$增大,$y$减小
则有最小正整数解$(x\%d_x+d_x)\%d_x$($d_x$为上文所提$db$)(若上式结果为0则最小正整数解为$d_x$),此时$y$取到两数均为正整数时的最大值
乘法逆元
$ab≡ba≡1(mod\ p)$
则称$b$是$mod\ p$意义下$a$的乘法逆元(定义了剩余系中的除法),记为$a^{-1}$,$a$的乘法逆元的意义是$a$在$mod\ p$剩余系下的倒数
乘法逆元的扩展欧几里得求法
给定$a,p$计算$b$使得$a*b≡1(mod\ p)$,等价于解方程$ax+py=1$($x,y$为整数变量,扩展欧几里得即可)
条件:$a,p$互质,即二元一次不定方程有解的条件,也即乘法逆元存在的充要条件
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
if(exgcd(a,mod,x,y)==1) //判断逆元是否存在 printf("%d\n",(x%mod+mod)%mod); //xj=即为a的逆元
时间复杂度$O(logn)$
费马小定理
假如$p$是质数,且$gcd(a,p)=1$,那么$a^{p-1}≡1(mod\ p)$
证明:
引理:集合$S1=\{1,2,3…p-1\},S2=\{1a,2a,3a…(p-1)a\}$,$S1$在模$p$意义下等于$S2$
故$(p-1)!≡a^{p-1}*(p-1)! (mod\ p)$,$(p-1)!$与$p$互素存在乘法逆元,故$a^{p-1}≡1(mod\ p)$。
乘法逆元的快速幂求法
若$p$是素数,由于$a^{p-1}≡1(mod\ p)$,故$a$的乘法逆元为$a^{p-2}$(就一个快速幂,所以这里就不再提供代码了)
条件:$a,p$互质且$p$为质数;时间复杂度:$O(logn)$
线性预处理乘法逆元
在$O(n)$的时间内求出$1..n$的逆元,条件:$1...n$均与$p$互素(为简化表达常表达为$p$为质数)
显然$1^{-1}≡1(mod\ p)$
设$p=\left\lfloor\dfrac{p}{k}\right\rfloor*k+r(0<r<k<p)$
可得$\left\lfloor\dfrac{p}{k}\right\rfloor*k+r≡0(mod\ p)$
两边同乘$k^{-1},r^{-1}$可得
$\left\lfloor\dfrac{p}{k}\right\rfloor*r^{-1}+k^{-1}≡0(mod\ p)$
可得$k^{-1}≡-\left\lfloor\dfrac{p}{k}\right\rfloor*(p\ mod\ k)^{-1}(mod\ p)$
由此可以实现递推求解逆元
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
a[1]=1; for(int i=2;i<=n;i++) { a[i]=-(p/i)*a[p%i]; a[i]=(a[i]%p+p)%p; }
积性函数
在数论中的积性函数:对于正整数$n$的一个算术函数$f(n)$,若$f(1)=1$,且当$a,b$互质时$f(ab)=f(a)f(b)$,在数论上就称它为积性函数。若对于某积性函数$f(n)$,就算$a,b$不互质,也有$f(ab)=f(a)f(b)$,则称它为完全积性的。
欧拉函数
对正整数$n$,欧拉函数$\varphi(n)$是小于等于$n$的正整数中与$n$互质的数的数目,例如$\varphi(8)=4(1,3,5,7)$
由欧拉函数的积性可得$φ(n)=\prodφ(p_i^{q_i})$
$φ(x=p^q)=p^{q-1}*(p-1)=x*(1-\frac{1}{p})$($p$为质数)
这个结论是显然的。
故$φ(x)=x*\prod\limits_{i=1}^n(1-\frac{1}{p_i})$
由$φ$的积性性质可由线性筛法$O(n)$计算$φ(1)...φ(n)$
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
p[1]=1,phi[1]=1; for(int i=2;i<=maxn-10;i++) { if(!p[i]) pr[++cnt]=i,phi[i]=i-1; for(int j=1;j<=cnt&&pr[j]*i<=maxn-10;j++) { p[i*pr[j]]=1; if(i%pr[j]==0) { phi[i*pr[j]]=phi[i]*pr[j]; break; } phi[i*pr[j]]=phi[i]*phi[pr[j]]; } }
欧拉定理
若$p,a$为正整数,且$p,a$互质,则:$a^{φ(p)}≡1(mod\ p)$
推论:$(a^{φ(p)})^{-1}=a^{φ(p)}$
推论:$a^b≡a^{b\%φ(p)}≡a^{b-k*φ(p)}(gcd(a,p)=1)$(等式两边同乘$(a^{φ(p)})^{-1}$)
扩展欧拉定理
$a^b≡\begin{cases}a^b,b<\varphi(p)\\a^{b\ mod\ \varphi(p)+\varphi(p)},b\geqslant\varphi(p)\end{cases}\pmod p$
不要求$a,p$互质
中国剩余定理(CRT)
中国剩余定理可用于求解同余方程组
$\begin{cases}x≡a_1(mod\ m_1)\\x≡a_2(mod\ m_2)\\.....\\x≡a_n(mod\ m_n)\end{cases}$
其中$m_1,m_2...m_n$为两两互质的整数
求解过程:
我们设$M=\prod\limits_{i=1}^nm_i\ ,\ M_i=\frac{M}{m_i}\ ,\ M_it_i≡1(mod\ m_i)\ (1<=i<=n)$
则该方程组必有一解$x_0=\sum\limits_{i=1}^na_iM_it_i$
可得通解$x=x_0+k*M$($k$为任意整数)
则有最小正整数解$x_{min}=(x\%M+M)\%M$,上式结果为0时$x_{min}=M$
证明:
$\begin{cases}a_jM_jt_j≡0(mod\ m_i)(i\ne j)\\a_jM_jt_j≡a_j(mod\ m_i)(i=j)\end{cases}$
故$\sum\limits_{i=1}^na_iM_it_i≡a_i(mod\ m_i)$
扩展中国剩余定理(EXCRT)
扩展中国剩余定理可用于求解同余方程组
$\begin{cases}x≡a_1(mod\ m_1)\\x≡a_2(mod\ m_2)\\.....\\x≡a_n(mod\ m_n)\end{cases}$
和中国剩余定理的区别在于,扩展中国剩余定理不要求$m_i$两两互质(很强的性质)
考虑递推求解,假设已经求出前$k-1$个方程组成的同余方程组的一个解$x$,设$M=lcm(m_1,m_2...m_{k-1})$
则前$k-1$的方程的通解为$x+i*M$($i$为任意整数)
那么对于加入第$k$个方程后的方程组
我们就是要求一个正整数$t$,使得$x+t*M≡a_k(mod\ m_k)$
即$t*M≡a_k-x(mod\ m_k)$
对于这个式子我们已经可以通过扩展欧几里得求解$t$
若该同余式无解,则整个方程组无解,若有,则前$k$个同余式组成的方程组的一个解为$x_k=x+t*M$
核心代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
scanf("%lld%lld%lld",&n,&a,&b); m=a,ans=b; for(int t,tmp,i=2;i<=n;ans%=m,i++) { scanf("%lld%lld",&a,&b); int g=exgcd(m,a,t,tmp),c=((b-ans)%a+a)%a; if(c%g) { printf("-1\n"); return 0; } t=mul(t,c/g,a/g),ans+=t*m,m=lcm(m,a); } ans=(ans%m+m)%m; printf("%lld\n",ans?ans:m);
大步小步算法(BSGS)
$BSGS$主要用于求解形如$a^x≡b(mod\ p)$的方程,其中$gcd(a,p)=1$
考虑最小解$x_0$,由欧拉定理的推论可得$a^x≡a^{x-k*φ(p)}$,则$x_0\in[0,φ(p)),k\in Z$
设$a^{x=ik-t}≡b(mod\ p)$,则必有$a^{ik}≡b*a^t(mod\ p)(0<t<=k,0<i<=\left\lceil\dfrac{φ(p)}{k}\right\rceil)$
枚举$t$,将$b*a^t\%p$存入哈希表(这里也可以使用$map$,不过会多带一个$log$),枚举$i$,在哈希表中查找即可得到第一个循环节中的所有解。第一个查找到的解即为最小解,对于第一个循环节中的任一解$x$,显然$x+k*φ(p),k\in Z^*$为该方程的解
时间复杂度$O(k+\dfrac{φ(p)}{k})$,由均值不等式可得$k=\sqrt{φ(p)}$时有最低时间复杂度$O(\sqrt{φ(p)})$
最小解求解代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
k=ceil(sqrt(phi)); int tp=1; for(int tmp=n*b%p,i=1;i<=k;i++,(tp*=b)%=p,(tmp*=b)%=p) mp[tmp]=i; for(int tmp=tp,i=1;i<=(phi+k-1)/k;i++,(tmp*=tp)%=p) if(mp[tmp]) { printf("%lld\n",i*k-mp[tmp]); return 0; } printf("no solution\n");
扩展BSGS
扩展$BSGS$主要用于求解形如$a^x≡b(mod\ p)$的方程
对于同余方程$a^x≡b(mod\ p)$,我们可以将其展开为$a*a^{x-1}+pk=b$的形式,可得$gcd(a,p)|b,x>0$为原方程有整数解的一个充分条件,这时讨论$x=0$等式是否成立,若成立则求得答案,否则若方程无解则同余方程无解
对于有解的方程,两边同除$gcd(a,p)$可得$\frac{a}{gcd(a,p)}*a^{x-1}+\frac{p}{gcd(a,p)}*k=\frac{b}{gcd(a,p)}$
转化为同余式可得$\frac{a}{gcd(a,p)}a^{x-1}≡\frac{b}{gcd(a,p)}(mod\ \frac{p}{gcd(a,p)})$
设$p'=\frac{p}{gcd(a,p)}$,若$gcd(a,p')=1$,按$BSGS$处理即可,否则重复执行上述过程
显然递归深度不超过$logp$,故可近似认为扩展$BSGS$与$BSGS$时间复杂度相同
完整代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<cmath> #include<map> #define int long long using namespace std; int c,a,p,b,k,ans,phi; map<int,int>mp; int gcd(int x,int y) { return y?gcd(y,x%y):x; } signed main() { while(1) { scanf("%lld%lld%lld",&a,&p,&b); if(!a) return 0; a%=p,mp.clear(),ans=0,c=1; int ps=p,bs=b; bool f=0; if(b>=p) { printf("No Solution\n"); continue; } if(!a) { printf("0\n"); continue; } for(int g=gcd(a,p);g!=1;g=gcd(a,p)) { if(c==b) { printf("%lld\n",ans),f=1; break; } if(b%g) { printf("No Solution\n"),f=1; break; } p/=g,(c*=a/g)%=p,b/=g,ans++; } if(f) continue; phi=p; for(int tmp=p,i=2;i*i<=tmp;i++) if(tmp%i==0) { while(tmp%i==0) tmp/=i; phi=phi*(i-1)/i; } k=ceil(sqrt(phi)); int tp=1; for(int tmp=b*a%p,i=1;i<=k;i++,(tp*=a)%=p,(tmp*=a)%=p) mp[tmp]=i; for(int tmp=tp,i=1;i<=(phi+k-1)/k;i++,(tmp*=tp)%=p) if(mp[tmp*c%p]) { printf("%lld\n",i*k-mp[tmp*c%p]+ans),f=1; break; } if(!f) printf("No Solution\n"); } return 0; }
阶和原根
定义:
阶:若$p,a$为正整数,且$p,a$互质,称满足同余式$a^x \equiv 1 \pmod p$的最小正整数$x$作$a$模$p$的阶,记作$ord_pa$
原根:若$p,a$为正整数,且$p,a$互质,且$ord_pa=\varphi(p)$,则称$a$为模$p$的原根
阶的性质:
1,$a^n≡1(mod\ p)$的充要条件为$\mathrm{ord}_pa|n$
推论:$ord_pa∣φ(p)$
2,对正整数$b$,若$a\equiv b\pmod p$,则$\text{ord}_pa=\text{ord}_pb$
3,设$m,n\in\mathbb N$,$a^m\equiv a^n\pmod p$的充要条件为 $m\equiv n\pmod {\text{ord}_pa}$
4,令$\text{ord}_pa=x$,则$1,a,a^2,\cdots,a^{x-1}$,模$p$两两不同余
5,对正整数$b$,若$ab\equiv1\pmod p$,则$\text{ord}_pa=\text{ord}_pb$
6,令$x=\text{ord}_pa$,对正整数$t$有$\text{ord}_pa^t=\frac{x}{\gcd(t,x)}$
7,对正整数$m$,若$m\mid p$,则$\text{ord}_ma\mid \text{ord}_pa$
8,对正整数$m$,若$m\bot p,a\bot m$,则 $\text{ord}_{mp}a=\text{lcm}(\text{ord}_{m}a,\text{ord}_{p}a)$
9,对正整数$b$,若$b\bot p$且$\mathrm{ord}_pa\bot \mathrm{ord}_pb$,则$\mathrm{ord}_pab=\mathrm{ord}_pa\times \mathrm{ord}_pb$
原根存在定理
原根存在的充要条件:模数为$2,4,p^k,2\times p^k$,其中$p$为奇素数,$k$为正整数
原根判定定理
$a$为模$m$的原根的充要条件是对于$\varphi(m)$的任意质因子$p$,必有 $a^{\frac{\varphi(m)}{p}}\not\equiv 1 \pmod m$,且$a^{\varphi(m)} \equiv 1 \pmod m$。
证明:由阶的性质1及其推论可得
原根的性质
1,设$a$为模$p$的原根,则集合$S=\{a^{s} \mid 1 \leq s \leq \varphi(p),s\bot\varphi(p)\}$给出模$p$的全部两两不同余的原根。模$p$的全部两两不同余的原根有 $\varphi(\varphi(p))$个。
2,最小原根是不大于$p^{\frac{1}{4}}$级别的
原根的求解
基于原根的性质,设计以下算法:依据原根判定定理枚举求解最小原根$a$;枚举指数$s$,若$s$与$\varphi(p)$互质,则$a^s\ mod \ p$为一个原根
完整代码
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<algorithm> #define int long long using namespace std; const int maxn=1e6+10; int t,n,phi[maxn],pr[maxn],g[maxn],cnt,ord[maxn]; bool p[maxn],f[maxn]; int gcd(int x,int y) { return y?gcd(y,x%y):x; } int quick_pow(int x,int y) { int ret=1; for(;y;y>>=1,(x*=x)%=n) if(y&1) (ret*=x)%=n; return ret; } bool check(int x) { for(int i=phi[n];i>1;i/=g[i]) if(quick_pow(x,phi[n]/g[i])==1) return 0; return 1; } signed main() { p[1]=1,phi[1]=1; for(int i=2;i<=maxn-10;i++) { if(!p[i]) pr[++cnt]=i,phi[i]=i-1,g[i]=i; for(int j=1;j<=cnt&&pr[j]*i<=maxn-10;j++) { p[i*pr[j]]=1,g[i*pr[j]]=pr[j]; if(i%pr[j]==0) { phi[i*pr[j]]=phi[i]*pr[j]; break; } phi[i*pr[j]]=phi[i]*phi[pr[j]]; } } f[2]=f[4]=1; for(int i=2;i<=cnt;i++) for(int j=pr[i];j<=maxn-10;j*=pr[i]) { f[j]=1; if(2*j<=maxn-10) f[2*j]=1; } scanf("%lld",&t); while(t--) { scanf("%lld",&n); if(!f[n]) { printf("0\n\n"); continue; } int a=1,tot=1; for(;!check(a)||quick_pow(a,phi[n])!=1;a++); ord[1]=a; for(int i=2;i<phi[n];i++) if(gcd(i,phi[n])==1) ord[++tot]=quick_pow(a,i); sort(ord+1,ord+tot+1); printf("%lld\n",tot); for(int i=1;i<=tot;i++) printf("%lld ",ord[i]); printf("\n"); } return 0; }
二次剩余
定义
二次剩余,俗称模意义开根,这里只介绍模数$p$为奇素数的解法
对于正整数$n$,若对于一个奇素数$p$,存在$x$满足$x^2\equiv n(mod \ p)$则称$n$是模$p$的二次剩余,否则称$n$是模$p$的非二次剩余
二次剩余的性质
1,对于模$p$的二次剩余$n$,该方程有两个不同解,且它们互为模意义下的相反数
2,任意一对相反数都对应一个二次剩余,而且这些二次剩余是两两不同的
3,二次剩余的数量恰为$\frac{p-1}{2}$,其他的非$0$数都是非二次剩余,数量也是$\frac{p-1}{2}$
4,二次剩余的逆元还是二次剩余,两个二次剩余的乘积还是二次剩余
欧拉准则
1,$n^{\frac{p-1}{2}} \equiv 1$与$n$是二次剩余等价
2,$n^{\frac{p-1}{2}} \equiv -1$与$n$是非二次剩余等价
$Cipolla$算法
首先找到一个$a$满足$a^2 - n$是非二次剩余,由于非二次剩余的数量是$\frac{p-1}{2}$,通过随机+根据欧拉准则检验的方式期望约$2$次可以找到这样一个$a$。
接下来定义$i^2 \equiv a^2 - n$。
但是$a^2 - n$不是二次剩余,怎么找得到这样一个$i$?
类比实数域到复数域的推广,定义这样一个$i$,将所有数表示为$A+Bi$的形式,
其中$A, B$都是模$p$意义下的数,类似于实部和虚部。
考虑证明$(a + i)^{p+1} \equiv n$
引理 1 :$i^p \equiv -i$
证明:$i^p \equiv i(i^2)^{\frac{p-1}{2}} \equiv i(a^2 - n)^{\frac{p-1}{2}} \equiv -i$
引理 2 :$(A + B)^p \equiv A^p + B^p$
证明:二项式定理展开后,由于$p$是质数,除了$C_p^0, C_p^p$外的组合数分子上的阶乘没法消掉,模$p$都会为$0$,剩下来的就是$C_p^0 A^0 B^p + C_p^p A^p B^0$
现在根据上述引理和费马小定理完成证明:
$(a + i)^{p+1} \equiv (a^p + i^p) (a + i) \equiv (a - i) (a + i) \equiv a^2 - i^2 \equiv n$
可以证明$(a + i)^{\frac{p+1}{2}}$的虚部一定为$0$,那么$(a + i)^{\frac{p+1}{2}}$即是一个解,其相反数是另一个解。
![](https://images.cnblogs.com/OutliningIndicators/ContractedBlock.gif)
#include<iostream> #include<cstdio> #include<cstdlib> #include<ctime> #define int long long using namespace std; struct c{ int x,y; }; int t,n,p,a,pw; int mod(int x) { return (x%p+p)%p; } c operator *(c &x,const c &y) { return (c){mod(x.x*y.x+pw*x.y%p*y.y),mod(x.x*y.y+y.x*x.y)}; } int quick_pow(c x,int y) { c ret;ret.x=1,ret.y=0; for(;y;y>>=1,x=x*x) if(y&1) ret=ret*x; return ret.x; } int pow(int x,int y) { int ret=1; for(;y;y>>=1,(x*=x)%=p) if(y&1) (ret*=x)%=p; return ret; } signed main() { srand(time(0)); scanf("%lld",&t); while(t--) { scanf("%lld%lld",&n,&p); if(pow(n,(p-1)/2)==p-1) { printf("Hola!\n"); continue; } else if(!n) { printf("0\n"); continue; } for(a=rand()%p;pow(pw=mod(a*a-n),(p-1)/2)!=p-1;a=rand()%p); int ans=quick_pow((c){a,1},(p+1)/2); if(p-ans<ans) ans=p-ans; printf("%lld %lld\n",ans,p-ans); } return 0; }