原根 (ex)BSGS 二次剩余 NTT
阶
由欧拉定理得 \(a^{\varphi(m)}\equiv1\space(\text{mod}\space m),(a,m)=1\) .
故满足 \(a^n\equiv1\space(\text{mod}\space m)\) 的最小 \(n\) 存在,称为 \(a\) 模 \(m\) 的阶,记作 \(\delta_m(a)\) .
-
\(a,a^2,...a^{\delta_m(a)}\) 模 \(m\) 两两不同余。
-
若 \(a^n\equiv 1\space(\text{mod}\space m)\),则 \(\delta_m(a)\mid n\) .
易证,否则与阶的最小性矛盾。
- 若 \((a,m)=(b,m)=1\),则 \(\delta_m(ab)=\delta_m(a)\delta_m(b)\leftrightarrow(\delta_m(a),\delta_m(b))=1\) .
必要性:
由于
有
故 \((\delta_m(a),\delta_m(b))=1\) .
充分性:
又因
故 \(\delta_m(ab)=\delta_m(a)\delta_m(b)\) .
- 若 \((a,m)=1\),则 \(\delta_m(a^k)=\frac{\delta_m(a)}{(\delta_m(a),k)}\) .
故 \(\delta_m(a^k)=\frac{\delta_m(a)}{(\delta_m(a),k)}\) .
原根
若 \((a,m)=1,\delta_m(a)=\varphi(m)\),则称 \(a\) 为模 \(m\) 的原根。
- 原根判定定理:若 \(m\ge3,(a,m)=1\),则 \(a\) 为模 \(m\) 的原根的充要条件是 \(\forall_{p\in\text{P},p\mid\varphi(m)},a^{\frac{\varphi(m)}{p}}\not\equiv1\space(\text{mod}\space m)\) .
必要显然,证充分。
反证法,若有 \(a^{\frac{\varphi(m)}{p}}\not\equiv1\space(\text{mod}\space m)\),假设存在一个 \(a\) 不是模 \(m\) 的原根。
则存在 \(t<\varphi(m),a^t\equiv1\space(\text{mod}\space m)\) .
由裴蜀定理得存在 \(k,x\) 满足 \(kt=x\varphi(m)+(t,\varphi(m))\) .
由欧拉定理得
矛盾,得证。
- 原根个数:\(\varphi(\varphi(m))\)(\(m\) 存在原根).
若有原根 \(g\),则
\(g^k\) 也是模 \(m\) 的原根即 \((k,\varphi(m))=1\) .
故原根有 \(\varphi(\varphi(m))\) 个。
- 原根存在定理:\(m\) 存在原根当且仅当 \(m=2,4,p^{\alpha},2p^{\alpha}\),其中 \(p\) 为奇素数。
这 part 太困难就不写了,直接记结论。
- 任意 \(m\) 的最小原根级别为 \(m^{0.25}\) .
求 \(n\) 的所有原根,从小到大记为 \(g_1,g_2,...g_c\) .
将所有下标为 \(d\) 为倍数的 \(g\) 输出。
\(T\le 10\),\(2\le n\le10^6\) .
判断一个数有没有原根,可以先预处理 \(p,p^{\alpha},2p^{\alpha}\),\(\text{O}(1)\) 判断即可。
若 \(g\) 为模 \(n\) 的原根:
预处理 \(\varphi(n)\) 的所有质因数,快速幂判断。
找到最小原根 \(g\),枚举其所有幂次 \(g^x,(x,\varphi(m))=1\) 即可。
#include<bits/stdc++.h>
#define ll long long
#define N 1000001
using namespace std;
int T,n,d,phi[N];
bool vis[N],flag[N];
vector<int>pr[N];
int gcd(int a,int b){
return !b?a:gcd(b,a%b);
}
void init(){
flag[2]=flag[4]=true;
for(int i=1;i<N;i++)phi[i]=i;
for(int i=2;i<N;i++){
if(vis[i])continue;
if(i!=2){
for(ll j=i;j<N;j*=i){
flag[j]=true;
if(j*2<N)flag[j*2]=true;
}
}
for(int j=i;j<N;j+=i){
vis[j]=true,phi[j]-=phi[j]/i;
pr[j].push_back(i);
}
}
}
int qpow(int k,int b,int p){
int ret=1;
while(b){
if(b&1)ret=1ll*ret*k%p;
k=1ll*k*k%p,b>>=1;
}
return ret;
}
int ans,g,rt[N],now,cur;
int main(){
scanf("%d",&T),init();
while(T--){
scanf("%d%d",&n,&d);
if(!flag[n]){
printf("0\n\n");
continue;
}
ans=phi[phi[n]],g=1,now=0;
printf("%d\n",ans);
if(n==2){
if(d==1)printf("1 ");
printf("\n");continue;
}
while(true){
if(gcd(g,n)!=1){
g++;continue;
}
bool fl=true;
for(int i=0,p;i<pr[phi[n]].size();i++){
p=pr[phi[n]][i];
if(qpow(g,phi[n]/p,n)==1){
fl=false;break;
}
}
if(fl)break;g++;
}
for(int i=1,cur=g;i<=phi[n];i++){
if(gcd(i,phi[n])==1)rt[++now]=cur;
cur=1ll*cur*g%n;
}
sort(rt+1,rt+1+ans);
for(int i=d;i<=ans;i+=d)
printf("%d ",rt[i]);
printf("\n");
}
return 0;
}
离散对数
对于模 \(n\) 的原根 \(g\) ,若 \(g^t\equiv a\space(\text{mod}\space n)\),称 \(a\) 在模 \(n\) 意义下,以 \(g\) 为底的对数为 \(t\) .
这个东西的 \(T=\varphi(n)\) .
多值函数关系如下:
主值为 \(\log_g a=t\) .
性质:
记号记为
换底公式
为使模 \(\varphi(n)\) 的除法能够进行,必须使用原根为底。
大步小步算法
BSGS,用于求解离散对数问题。
求 \(a^x\equiv b\space(\text{mod}\space p)\) 的解,其中 \(a\perp p\) .
令 \(x=A\lceil\sqrt{p}\rceil-B\space(0\le A,B\le\lceil\sqrt{p}\rceil)\),则 \(a^{A\lceil\sqrt{p}\rceil}\equiv ba^B\space(\text{mod}\space p)\) .
枚举所有 \(ba^B\),用 \(\text{map}\) 存,计算 \(a^{A\lceil\sqrt{p}\rceil}\),寻找与之相符的数,可以得到所有的 \(x=A\lceil\sqrt{p}\rceil-B\) .
复杂度 \(\text{O}(\sqrt{V}\log\sqrt{V})\) .
P3846 [TJOI2007] 可爱的质数/【模板】BSGS
求 \(a^x\equiv b\space(\text{mod}\space p)\) 的最小正整数解,无解输出 no solution
.
这里 \(p\in\text{P}\),值域 \(2^{31}\) .
#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll qpow(ll k,ll b,ll p){
ll ret=1;k%=p;
while(b){
if(b&1)(ret*=k)%=p;
(k*=k)%=p,b>>=1;
}
return ret;
}
ll a,b,p;
ll BSGS(ll a,ll b,ll p){
map<ll,ll>hs;
hs.clear(),b%=p;
ll B=sqrt(p)+1;
for(ll i=0,cur=b;i<B;i++)
hs[cur]=i,(cur*=a)%=p;
a=qpow(a,B,p);
if(!a)return !b?1:-1;
for(ll i=1,cur=1;i<=B;i++){
(cur*=a)%=p;
ll j=hs.find(cur)==hs.end()?-1:hs[cur];
if(j>=0&&i*B-j>=0)return i*B-j;
}
return -1;
}
int main(){
cin>>p>>a>>b;
ll ans=BSGS(a,b,p);
if(ans==-1)printf("no solution\n");
else printf("%lld\n",ans);
return 0;
}
二次剩余
若 \(a\) 非 \(p\) 的倍数且模 \(p\) 同余于某个数的平方,则 \(a\) 为模 \(p\) 的 二次剩余 ,否则 \(a\) 为模 \(p\) 的 二次非剩余 。
对二次剩余求解:
即求模意义下的开平方运算。只对 \(p\) 为 奇素数 的情况讨论。
\(\text{Legendre}\) 符号
记为
\(p\mid a\Rightarrow (\frac{a}{p})=0\) .
\(a\) 为模 \(p\) 的二次剩余 \(\Rightarrow(\frac{a}{p})=1\) .
\(a\) 为模 \(p\) 的二次非剩余 \(\Rightarrow(\frac{a}{p})=-1\) .
\(\text{Euler}\) 判别准则
若 \(p\not|\space a\),则
二次剩余可解同余 \(1\),否则 \(-1\) .
引理:令 \(a\equiv g^k\) ,那么二次剩余有解当且仅当 \(k\) 为偶数。
充分:若 \(x^2\equiv a\space(\text{mod}\space p)\) 对于某个 \(g^l\) 有解,则
\(2\mid p-1\Rightarrow 2\mid k\) .
必要:假设 \(k\) 为偶数,则
由于 \(g^{p-1}\equiv 1\space(\text{mod}\space p)\) ,且 \(\forall 1\le k<p-1,g^k\not\equiv1\space(\text{mod}\space p)\),有
那么
所以 \(a^{(p-1)/2}\equiv1\space(\text{mod}\space p)\) 时有解。
由引理得无解时 \(k\) 为奇数,则
得到 \(\text{Euler}\) 判别准则,且 \(\text{Legendre}\) 符号为完全积性函数。
二次剩余和二次非剩余的数量
感性方法。
对于两个不同的解 \(x_0,x_1\),可得 \((x_0-x_1)(x_0+x_1)\equiv 0\space(\text{mod}\space p)\) .
故这两个解为相反数,它们对应一组二次剩余,故 二次剩余与二次非剩余的数量都是 \(\frac{p-1}{2}\) .
如上,求解
\(\text{Cipolla}\) 算法
找到 \(a\) 满足 \(a^2-n\) 为二次非剩余,定义 \(i^2=a^2-n\),即将所有数在类似复数域上表示为 \(A+Bi\) 的形式。
证明:\((a+i)^{p+1}\equiv n\space(\text{mod}\space p)\) .
\(\text{Lemma 1}\):\(i^p\equiv -i\) .
得证。
\(\text{Lemma 2}\):\((A+B)^p\equiv A^p+B^p\) .
由二项式定理可以消去非 \(0\) 和 \(p\) 项上的系数,剩下的就是
于是
那么 \((a+i)^{\frac{p+1}{2}}\) 即为一个解。
- \((a+i)^{\frac{p+1}{2}}\) 的虚部为 \(0\) .
反证,假设存在 \((A+Bi)^2\equiv n\) 且 \(B\not=0\),那么
左式虚部为 \(0\),故右式的虚部也为 \(0\),即 \(AB\equiv 0\) .
因为 \(B\not=0\) 所以 \(A\equiv 0\),有 \((Bi)^2\equiv n\) .
得到 \(i^2\equiv nB^{-2}\),故 \(B^2\) 和 \(B^{-2}\) 都是二次剩余,乘上二次剩余 \(n\) 后仍为二次剩余,与 \(i^2\) 为二次非剩余矛盾。
- 如何找二次非剩余?
在 \(\lbrack0,p)\) 中随机一个 \(a\),用欧拉准则判定 \(a\) 是否为非二次剩余即可,期望 \(2\) 次找到。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
ll read(){
ll x=0,w=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
while(isdigit(ch))x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x*w;
}
ll T,w,n,p;
struct num{
ll x,y;
num operator*(const num &a)const{
num ret;ret.x=ret.y=0;
ret.x=(x*a.x%p+y*a.y%p*w%p)%p;
ret.y=(x*a.y%p+y*a.x%p)%p;
return ret;
}
};
ll qpow(ll k,ll b,ll p){
ll ret=1;
while(b){
if(b&1)(ret*=k)%=p;
(k*=k)%=p,b>>=1;
}
return ret;
}
ll poww(num k,ll b,ll p){
num ret=(num){1,0};
while(b){
if(b&1)ret=ret*k;
k=k*k,b>>=1;
}
return ret.x;
}
ll Cipolla(ll n,ll p){
n%=p;ll a;
if(qpow(n,(p-1)/2,p)==p-1)return -1;
while(true){
a=rand()%p,w=(a*a%p-n+p)%p;
if(qpow(w,(p-1)/2,p)==p-1)break;
}
num x=(num){a,1};
return poww(x,(p+1)/2,p);
}
int main(){
srand(time(0)),T=read();
while(T--){
n=read(),p=read();
if(!n){
printf("0\n");
continue;
}
ll ans1=Cipolla(n,p),ans2=p-ans1;
if(ans1==-1){
printf("Hola!\n");
continue;
}
if(ans1>ans2)swap(ans1,ans2);
if(ans1==ans2)printf("%lld\n",ans1);
else printf("%lld %lld\n",ans1,ans2);
}
return 0;
}
求 \(a^x\equiv b\space(\text{mod}\space p)\) 的最小非负整数解。
\(a,p\) 不一定互质 。
设 \(d_1=(a,p)\),若 \(d_1\not|\space b\) 显然无解。
将方程除以 \(d_1\) 得到
若 \(a\) 与 \(\frac{p}{d_1}\) 仍不互质,继续进行上述操作。
记 \(D=\prod_{i=1}^{k}d_i\),方程变为
\(a\perp\frac{p}{D}\rightarrow a^k\perp\frac{p}{D}\) ,这是一个普通的 BSGS 问题,将答案加上 \(k\) 即可。
对于所有 \(0\le i\le k\),暴力判断原式是否成立。
细节多点看一下。
#include<bits/stdc++.h>
#define ll long long
using namespace std;
int read(){
int x=0,w=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
while(isdigit(ch))x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x*w;
}
ll gcd(ll a,ll b){
return !b?a:gcd(b,a%b);
}
ll phi(ll k){
ll ret=k;
for(ll i=2;i*i<=k;i++)
if(k%i==0){
ret-=ret/i;
while(k%i==0)k/=i;
}
if(k>1)ret-=ret/k;
return ret;
}
ll qpow(ll k,ll b,ll p){
ll ret=1;k%=p;
while(b){
if(b&1)(ret*=k)%=p;
(k*=k)%=p,b>>=1;
}
return ret;
}
ll T,a,b,p;
ll pw,k,tp;
map<ll,ll>hs;
ll BSGS(ll a,ll b,ll p,ll pw){
hs.clear();
ll B=sqrt(p)+1;
for(ll i=0,cur=b;i<B;i++)
hs[cur]=i,(cur*=a)%=p;
a=qpow(a,B,p);
if(!a)return !b?1:-1;
for(ll i=1,cur=pw;i<=B;i++){
(cur*=a)%=p;
ll j=hs.find(cur)==hs.end()?-1:hs[cur];
if(j>=0&&i*B-j>=0)return i*B-j;
}
return -1;
}
ll exBSGS(ll a,ll b,ll p){
a%=p,b%=p;
if(b==1||p==1)return 0;
if(b%gcd(a,p))return -1;
k=0,pw=1;
while((tp=gcd(a,p))!=1){
if(b%tp)return -1;
k++,p/=tp,b/=tp;
(pw*=a/tp)%=p;
if(pw==b)return k;
}
ll ans=BSGS(a,b,p,pw);
return ans==-1?-1:ans+k;
}
int main(){
while(true){
a=read(),p=read(),b=read();
if(!a&&!p&&!b)break;
ll ans=exBSGS(a,b,p);
if(ans==-1)printf("No Solution\n");
else printf("%lld\n",ans);
}
return 0;
}
NTT
叫做(快速)数论变换,理解为在环上的 \(\text{FFT}\) .
对于质数 \(p=qn+1,n=2^k\),其原根 \(g\) 满足 \(g^{qn}\equiv1\space(\text{mod}\space p)\) .
一般有
将 \(g_n=g^q\space(\text{mod}\space p)\) 视作 \(\omega_n\) 的等价,可以发现
-
\(\omega_n^0=\omega_n^n\equiv1,\omega_n^x\equiv(\omega_n^1)^x\space(\text{mod}\space p)\) .
-
\(\omega_{2n}^{2x}\equiv\omega_n^x\space(\text{mod}\space p)\) .
-
\(\omega_{n}^{x}\equiv-\omega_n^{x+\frac{n}{2}}\space(\text{mod}\space p)\) .
证明:两边同除 \(\omega_n^x\) 可得 \(1\equiv-\omega_{n}^{\frac{n}{2}}\),由 \(\omega_{n}^{\frac{n}{2}}\equiv-1\space(\text{mod}\space p)\) 可证原式。
然后进行分治:
设
故 \(F(x)=A(x^2)+x\cdot B(x^2)\)
代入 \(\omega_{n}^{x}\) 有
代入 \(\omega_{n}^{x+\frac{n}{2}}\) 有
发现只需扫过 \(\omega_{n}^{0}\) 至 \(\omega_{n}^{\frac{n}{2}-1}\) 就可以得到另一半。
所以分治复杂度 \(\text{O}(n\log n)\) .
- \(A(x),B(x)\) 按照下标二进制最后一位划分,最后递归到的位置是翻转后的二进制,可以预处理。
没太搞懂贺下板子吧。
P1919 【模板】A*B Problem 升级版(FFT 快速傅里叶变换)
大数乘法,位数 \(10^6\) .
取模数 \(P=998244353\),板子背了就好。
#include<bits/stdc++.h>
#define ll long long
#define N 1<<21
#define P 998244353
using namespace std;
int read(){
int x=0,w=1;char ch=getchar();
while(!isdigit(ch)){if(ch=='-')w=-1;ch=getchar();}
while(isdigit(ch))x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
return x*w;
}
int qpow(int k,int b){
int ret=1;
while(b){
if(b&1)ret=1ll*ret*k%P;
k=1ll*k*k%P,b>>=1;
}
return ret;
}
int A[N],B[N],C[N];
char a[N],b[N];
int n,lim=1,r[N];
int m,k,gn,tp,inv;
void ntt(int *x,int lim,int opt){
for(int i=0;i<lim;i++)
if(r[i]<i)swap(x[i],x[r[i]]);
for(int m=2;m<=lim;m<<=1){
k=m>>1,gn=qpow(3,(P-1)/m);
for(int i=0;i<lim;i+=m){
for(int j=0,g=1;j<k;j++,g=1ll*g*gn%P){
tp=1ll*x[i+j+k]*g%P;
x[i+j+k]=(x[i+j]-tp+P)%P;
x[i+j]=(x[i+j]+tp)%P;
}
}
}
if(opt==-1){
reverse(x+1,x+lim),inv=qpow(lim,P-2);
for(int i=0;i<lim;i++)
x[i]=1ll*x[i]*inv%P;
}
}
int main(){
scanf("%s",&a),n=strlen(a);
for(int i=0;i<n;i++)A[i]=a[n-i-1]-'0';
while(lim<(n<<1))lim<<=1;
scanf("%s",&b),n=strlen(b);
for(int i=0;i<n;i++)B[i]=b[n-i-1]-'0';
while(lim<(n<<1))lim<<=1;
for(int i=0;i<lim;i++)
r[i]=(r[i>>1]>>1)+(i&1)*(lim>>1);
ntt(A,lim,1),ntt(B,lim,1);
for(int i=0;i<lim;i++)
C[i]=1ll*A[i]*B[i]%P;
ntt(C,lim,-1);
int len=0;
for(int i=0;i<lim;i++){
if(C[i]>=10){
len=i+1;
C[i+1]+=C[i]/10,C[i]%=10;
}
if(C[i])len=max(len,i);
}
while(C[len]>=10)
C[len+1]=C[len]/10,C[len]%=10,len++;
for(int i=len;i>=0;i--)
putchar(C[i]+'0');
printf("\n");
return 0;
}