初等数论（Ⅲ）：高次同余阶和原根相关

Posted 2023-05-27 bloodstalk

高次同余

前言

关于高次同余方程，有 \\(a^x \\equiv b(\\textmod \\ p)\\) 和 \\(x^a \\equiv b(\\textmod \\ p)\\) 两种类型，后者计算起来较为麻烦，下文就分别记述这两种高次同余方程。

离散对数问题

离散对数问题是在模 \\(p\\) 意义下求解 \\(\\log_ab\\)，这等价于形如

\\[a^x \\equiv b(\\textmod \\ p) \\]

的高次同余方程，其中 \\(x\\) 即为 \\(\\log_ab\\)，\\(x\\) 是一个非负整数。

当 \\(a \\perp b\\) 时，我们可以采用 BSGS 算法解决；当 \\(a \\not\\perp b\\) 时，可以采用 exBSGS 算法求解。

BSGS

大步小步算法，英文名为 Baby Step Giant Step，简称 BSGS。适用于 \\(a\\perp p\\) 的情况。

算法流程

由扩展欧拉定理得

\\[a^x \\equiv a^x \\mod \\varphi(p)(\\textmod \\ p) \\]

所以 \\(a^x\\) 模 \\(p\\) 意义下的循环节就是 \\(\\varphi(p)\\)，又因 \\(\\varphi(p) < p\\)，所以考虑 \\(x\\in[0,p]\\) ，一定能找到最小整数 \\(x\\)。自此，问题缩小到了 \\(O(p)\\) 级别。如果 \\(p\\leq 2^32\\)，那就超时了。这时候就要请出我们今天的主角----根号平衡闪亮登场！

我们令 \\(x = im-j\\)，其中 \\(m=\\lceil \\sqrtp \\ \\rceil , i\\in[1,m],j\\in[0,m-1]\\)。这样分配，\\(im-j\\) 就能取遍 \\([1,p]\\) 了，如果 \\(x=0\\)，说明 \\(b=1\\)，特判一下就行了。

我们把上式转换一下，则有

\\[a^im-j \\equiv b(\\textmod \\ p) \\]

因为 \\(a\\perp p\\) ，所以

\\[(a^m)^i \\equiv ba^j(\\textmod \\ p) \\]

是一个等价转换（除回去模数不变）。

然后枚举 \\(i,j\\)。

1. 先枚举 \\(j\\)，把 \\((ba^j \\ \\textmod \\ p,j)\\) 插入一个 Hash 表。如果 \\(ba^j \\ \\textmod \\ p\\) 出现相等的情况，为了求最小解，显然用更大的 \\(j\\) 替代小的解。

2. 然后枚举 \\(i\\) ，计算 \\((a^m)^i \\ \\textmod \\ p\\)，到 Hash 表中判断是否有相等的 key，因为 \\(i\\times m\\) 的变化量是比 \\(j\\) 大的，所以我们找到第一个就可以结束程序，最小的 \\(x=im-j\\) 就出来了。

这样，枚举 \\(i,j\\) 的次数都是 \\(\\sqrt p\\) 的，总算法的复杂度也就是 \\(O(\\sqrtp\\ )\\) 的。

map <int,int> Hash;

signed main()

	mod = read() , a = read(), b = read();
	if(b == 1)  return printf("0"),0; //特判一下
	m = ceil(sqrt(mod));
	Hash[b] = 0/*j取0,ba^j就是b*/ , t = b;
	for(re int i=1;i<m;i++)//枚举ba^j
	
		t = t * a % mod;
		Hash[t] = i;
	
    
    
	val = ksm(a,m) , t = 1;
	for(re int i=1;i<=m;i++)//枚举(a^m)^i
	
		t  = t * val % mod;
		if(Hash.count(t))   return printf("%lld",i*m-Hash[t]),0; //有值直接输出
	
	printf("no solution");
	return 0;

exBSGS

当 \\(a\\not\\perp p\\) 的时候，我们的 exBSGS 就要出场了。

化未知为已知，我们思考怎么操作能让 \\(a,p\\) 互质。

首先，原方程可以写作

\\[a \\cdot a^x-1 \\equiv b(\\textmod \\ p) \\]

的形式，这就是把 \\(a\\) 提出来的一个过程。这也等价于求 \\(a\\cdot a^x-1 + py = b\\) 的解。

令 \\(d_1 = \\gcd(a,p)\\)。如果 \\(d_1\\nmid b\\)，则原方程无解。

否则由同余的性质得，同余方程就变成

\\[\\fracad_1a^x-1 \\equiv \\fracbd_1(\\textmod \\ \\fracpd_1) \\]

这启发我们要不断提取 \\(a\\) 出来，直到 \\(a\\) 和模数互质。

设 \\(d_2 = \\gcd(a,\\fracpd_1)\\)，然后重复上述操作，直到互质。

当 \\(a\\perp \\fracpd_1\\dots d_k\\) 时，我们就可以套用 BSGS 了。设 \\(D = \\prod_i=1^kd_i\\)，原方程就变成了形如

\\[\\fraca^kDa^x-k \\equiv \\fracbD(\\textmod \\ \\fracpD) \\]

的形式。因为 \\(a\\perp \\fracpD\\)，所以 \\(\\fraca^kD\\perp \\fracpD\\)，那么就可以求个逆元，把 \\(\\fraca^kD\\) 消掉，然后再套一个 BSGS 就可以了。

注意，BSGS 结束后算出来的是 \\(x-k\\) 而非 \\(x\\)，别忘了把 \\(k\\) 加上。

实际上，\\(\\fraca^kD\\) 也可以不用求逆元，只要在 BSGS 枚举 \\(i\\) 的时候把初值设为 \\(\\fraca^kD\\) 即可。

code

il int BSGS(int a,int b,int mod,int k)

	Hash.clear();
	m = ceil(sqrt(mod));
	Hash[b] = 0 , t = b;
	for(re int i=1;i<m;i++)
	
		t = 1ll * t * a % mod;
		Hash[t] = i;
	
	val = ksm(a,m,mod) , t = k;//初值设为k
	for(re int i=1;i<=m;i++)
	
		t  = 1ll * t * val % mod;
		if(Hash.count(t))   return i*m-Hash[t]; 
	
	return -1;


il void exBSGS()

	a = read() , mod = read() , b = read();
	if(!a && !mod && !b) exit(0);
	a %= mod , b %= mod;//先模再特判
	if(b == 1 || mod == 1)  puts("0"); return ; 
	cnt = d = 0 , k = 1;
	while((d=gcd(a,mod)) != 1)
	
		if(b % d)  puts("No Solution"); return; //无解情况
		cnt++;//不知道为什么这个放在if上面就wa一个点，很怪
		b /= d , mod /= d;//逐渐往下找
		k = 1ll * k * (a / d) % mod;//a^x-k前的系数
		if(k == b)  printf("%d\\n",cnt); return ; //k==b说明a^x-k = 1,x=k
	
	ans = BSGS(a,b,mod,k);
	if(ans == -1) puts("No Solution");
	else printf("%d\\n",ans+cnt);

解高次同余方程51nod1038 X^A Mod P

1038 X^A Mod P 

基准时间限制：1 秒 空间限制：131072 KB 分值: 320

X^A mod P = B，其中P为质数。给出P和A B，求< P的所有X。

例如：P = 11，A = 3，B = 5。

3^3 Mod 11 = 5

所有数据中，解的数量不超过Sqrt(P)。

 

Input

第1行：一个数T，表示后面用作输入测试的数的数量。（1 <= T <= 100)
第2 - T + 1行：每行3个数P A B，中间用空格隔开。(1 <= A, B < P <= 10^9, P为质数)

Output

共T行，每行包括符合条件的X，且0 <= X < P，如果有多个，按照升序排列，中间用空格隔开。如果没有符合条件的X，输出：No Solution。所有数据中，解的数量不超过Sqrt(P)。

Input示例

3
11 3 5
13 3 1
13 2 2

Output示例

3
1 3 9
No Solution

题解

解高次同余方程时，先求出来原根

再用原根代换方程两边（要用到BSGS）

转化成一次剩余问题用EX_GCD求解即可

P.S. BSGS用map会T，所以我恬不知耻地贴了一份

代码

//by 减维
#include<set>
#include<map>
#include<queue>
#include<ctime>
#include<cmath>
#include<bitset>
#include<vector>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define ll long long
#define il inline
#define rg register
#define db double
#define mpr make_pair
#define maxn 100005
#define inf (1<<30)
#define eps 1e-8
#define pi 3.1415926535897932384626L
using namespace std;

inline int read()
{
    int ret=0;bool fla=0;char ch=getchar();
    while((ch<\'0\'||ch>\'9\')&&ch!=\'-\')ch=getchar();
    if(ch==\'-\'){fla=1;ch=getchar();}
    while(ch>=\'0\'&&ch<=\'9\'){ret=ret*10+ch-\'0\';ch=getchar();}
    return fla?-ret:ret;
}

int t,num,cnt;
ll a,p,b,g,phip,pri[maxn],ans[maxn];
vector<ll> phipri;
map<ll,ll> mp;
bool pd[maxn];

il ll ksm(ll x,ll y,ll mod)
{
    ll ret=1;x%=mod;
    for(;y;y>>=1,x=x*x%mod)
        if(y&1) ret=ret*x%mod;
    return ret;
}

ll exgcd(ll a,ll b,ll &x,ll &y)
{
    if(b==0){x=1;y=0;return a;}
    ll gcd=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return gcd;
}

il void pre()
{
    for(int i=2;i<=maxn-5;++i)
    {
        if(!pd[i]) pri[++num]=i;
        for(int j=1;j<=num&&i*pri[j]<=maxn-5;++j)
        {
            pd[i*pri[j]]=1;
            if(i%pri[j]==0) break;
        }
    }
}

il ll getphi(ll x)
{
    ll ret=x;
    for(int i=1;pri[i]*pri[i]<=x;++i)
        if(x%pri[i]==0)
        {
            ret=ret/pri[i]*(pri[i]-1);
            while(x%pri[i]==0) x/=pri[i];
        }
    if(x>1) ret=ret/x*(x-1);
    return ret;
}

il bool check(ll g,ll x,ll p)
{
    int siz=phipri.size();
    for(int i=0;i<siz;++i)
        if(ksm(g,x/phipri[i],p)==1) return 0;
    return 1;
}

il ll getg(ll x,ll p)
{
    ll tmp=x;
    phipri.clear();
    for(int i=2;i*i<=x;++i)
        if(x%i==0)
        {
            phipri.push_back(i);//printf("%lld ",pri[i]);
            while(x%i==0) x/=i;
        }
    if(x!=1) phipri.push_back(x);//,printf("%lld ",x);
    //puts("");
    x=tmp;
    ll gen=1;
    while(1)
    {
        if(check(gen,x,p)) return gen;
        gen++;
    }
}

struct sa{
    long long x;
    int id;
    bool operator<(const sa &b)const{
        if (x == b.x) return id < b.id;
        return x<b.x;
    }
}rec[100500];
//用rec存离散对数
long long bsgs(long long x,long long n,long long m){
    int s=(int)(sqrt((double)m+0.5));
    while((long long)s*s<=m)s++;
    long long cur=1;
    sa tmp;
    for(int i=0;i<s;i++){
        tmp.x=cur,tmp.id=i;
        rec[i]=tmp;
        cur=cur*x%m;
    }
    sort(rec,rec+s);
    //这里不能用map查找比较慢，采用排序二分就快了
    long long mul= ksm(cur, m - 2, m) % m;
    //这里有的方法是在下面的循环里求解快速幂，但本题是不行的  要在循环外面弄，保证时间
    cur=1;

    for(long long i=0;i<s;i++){
        long long more=n*cur%m;
        tmp.x=more,tmp.id=-1;
        int j=lower_bound(rec,rec+s,tmp)-rec;
        if(rec[j].x==more){
            return i*s+rec[j].id;
        }
        cur=cur*mul%m;
    }
    return -1;
}

int main()
{
    t=read();
    pre();
    while(t--){
        p=read(),a=read(),b=read();
        phip=p-1;cnt=0;
        g=getg(phip,p);
        //printf("%lld ",g);
        ll c=bsgs(g,b,p);
        //printf("%lld ",c);
        if(c==-1){puts("No Solution");continue;}
        ll x,y;
        ll gcd=exgcd(a,phip,x,y);
        //printf("%lld ",gcd);
        if(c%gcd!=0){puts("No Solution");continue;}
        x=x*(c/gcd)%phip;
        ll delt=phip/gcd;
        for(int i=0;i<gcd;++i)
        {
            x=((x+delt)%phip+phip)%phip;
            ans[++cnt]=ksm(g,x,p);
        }
        sort(ans+1,ans+cnt+1);
        for(int i=1;i<=cnt;++i) printf("%lld ",ans[i]);puts("");
    }
    return 0;
}