<更新提示>
<第一次更新>
<正文>
高次同余方程
一般来说,高次同余方程分(a^x equiv b(mod p))和(x^a equiv b(mod p))两种,其中后者的难度较大,本片博客仅将介绍第一类方程的解决方法。
给定(a,b,p),其中(gcd(a,p)=1),求方程(a^x equiv b(mod p))的最小非负整数解。
普通分析和朴素算法
先介绍一下欧拉定理:
如果正整数(a),(p)互质,则(a^{phi(p)}equiv1(mod p))。
注意到题中所给的条件符合欧拉定理的前提条件,所以我们得知(a^{phi(p)}equiv 1(mod p)),又因为(phi(p)leq p-1),(a^0equiv 1(mod p)),所以在(0)到(p-1)范围内,模(p)意义下(a)的幂一定形成了至少一次循环。也就是说,我们只要在(0)到(p-2)范围内枚举解(x),就一定能找到解,用快速幂判断即可求解该方程,时间复杂度(O(plog_2p))。
Baby Step Giant Step算法
对于求解该类高次同余方程,(Baby Step Giant Step)可以优化枚举算法,在(O(sqrt p))的时间内求解该方程。
我们设解(x=i*t-j),其中(0leq j leq t-1),则方程为(a^{i*t+j}equiv b(mod p)⇔(a^t)^iequiv b*a^j(mod q))。
对于所有的(jin[0,t-1]),将(b*a^j mod p)的值插入(hash)表。
枚举(i)的所有可能取值([0,frac{p}{t}]),计算出((a^t)^i mod p)的值,并在(hash)表中检索该值,如果该值存在,就可以用(i*t-j)更新答案。
该算法可以保证枚举到(x)在([0,p-2])范围内的每一种情况,其时间复杂度为(O(frac{p}{t}+p)),该函数为经典的对勾函数,当参数(t)取(sqrt p)时,其函数值最小,保证其时间复杂度为(O(sqrt p))。
以下给出求解方程最小非负整数解的程序,无解时返回(-1)。
(Code:)
inline long long Baby_Step_Giant_Step(long long a,long long b,long long p)
{
map < long long , long long > hash;hash.clear();
b%=p;long long t=(long long)sqrt(p)+1,mul=1;
for (int j=0;j<t;j++)
{
hash[ mul * b % p ] = j;
mul = mul * a % p;
}
if (a%p==0)return b==0 ? 1 : -1;
a = 1;
for (int i=0;i<=t;i++)
{
long long j = ( hash.find(a) != hash.end() ? hash[a] : -1 );
if (j>=0&&i*t-j>=0)return i*t-j;
a = a * mul % p;
}
return -1;
}
计算器
Description
你被要求设计一个计算器完成以下三项任务:
1、给定y,z,p,计算Y^Z Mod P 的值;
2、给定y,z,p,计算满足xy≡ Z ( mod P )的最小非负整数;
3、给定y,z,p,计算满足Y^x ≡ Z ( mod P)的最小非负整数。
Input Format
输入包含多组数据。
第一行包含两个正整数T,K分别表示数据组数和询问类型(对于一个测试点内的所有数据,询问类型相同)。
以下行每行包含三个正整数y,z,p,描述一个询问。
Output Format
对于每个询问,输出一行答案。对于询问类型2和3,如果不存在满足条件的,则输出“Orz, I cannot find x!”,注意逗号与“I”之间有一个空格。
Sample Input
3 1
2 1 3
2 2 3
2 3 3
Sample Output
2
1
2
解析
同余方程模板题。对于询问(1),快速幂解决,对于询问(2),扩展欧几里得算法解决,对于询问(3),(BSGS)算法解决。
对于第三个任务,一个细节要注意的是,题中保证了(p)为质数,却没保证(a)不是(p)的倍数。所以要特判(a)为(p)倍数的情况。
(Code:)
#include<bits/stdc++.h>
using namespace std;
#define mset(name,val) memset(name,val,sizeof name)
long long t,k;
inline void input(void)
{
scanf("%lld%lld",&t,&k);
}
inline long long power(long long a,long long b,long long p)
{
long long res=1;
while (b)
{
if (1&b)res*=a,res%=p;
b>>=1;
a*=a,a%=p;
}
return res;
}
inline long long Exeuclid(long long a,long long &x,long long b,long long &y,long long c)
{
if (b==0){x=c/a,y=0;return a;}
else
{
long long p=Exeuclid(b,x,a%b,y,c);
long long x_=x,y_=y;
x=y_;y=x_-a/b*y_;
return p;
}
}
inline long long Euclid(long long a,long long b)
{
return b ? Euclid(b,a%b) : a;
}
inline long long Baby_Step_Giant_Step(long long a,long long b,long long p)
{
map < long long , long long > hash;hash.clear();
b%=p;long long t=(long long)sqrt(p)+1,mul=1;
for (int j=0;j<t;j++)
{
hash[ mul * b % p ] = j;
mul = mul * a % p;
}
if (a%p==0)return b==0 ? 1 : -1;
a = 1;
for (int i=0;i<=t;i++)
{
long long j = ( hash.find(a) != hash.end() ? hash[a] : -1 );
if (j>=0&&i*t-j>=0)return i*t-j;
a = a * mul % p;
}
return -1;
}
inline void solve(void)
{
long long y,z,p;
for (int i=1;i<=t;i++)
{
scanf("%lld%lld%lld",&y,&z,&p);
if (k==1)
printf("%lld
",power(y,z,p));
if (k==2)
{
if (z%Euclid(y,p))
{
printf("Orz, I cannot find x!
");
continue;
}
long long x_,y_;
Exeuclid(y,x_,p,y_,z);
printf("%lld
",(x_%p+p)%p);
}
if (k==3)
{
long long ans=Baby_Step_Giant_Step(y,z,p);
if (ans==-1)
printf("Orz, I cannot find x!
");
else printf("%lld
",ans);
}
}
}
int main(void)
{
input();
solve();
return 0;
}
ExBSGS算法
顾名思义,该算法就是(Baby Step Giant Step)的拓展。对于方程(a^xequiv b(mod p)),当(a,p)不互质的时候,我们可以使用该算法进行求解。
不难发现,当互质条件失去后,无解的情况也随之增加。经过严格的数学证明可以得到:当(gcd(a,p) ot|b)且(b ot=1)时,该方程无自然数解(转换成扩欧方程的有解性判定,由裴蜀定理可得)。
(b=1)的情况是需要特判的,(x)最小非负整数解为(0)。除了无解和该情况以外,我们可以用如下方法求解:
令(d=gcd(a,p)),将原式转换为 :$$a^{x-1} frac{a}{d} equiv frac{b}{d} (mod frac{p}{d})$$
设(x_1=x-1,k_1=frac{a}{d},b_1=frac{b}{d},p_1=frac{p}{d}),则
若能够求解该方程,则原方程的解(x=x_1+1)。
记录下(k),则这个方程的结构是可以进行同样的操作的。如此,不断缩小模数(p),直到(gcd(a,p)=1),得到方程$$k_1k_2...k_na^{x_n}equiv b_n(mod p_n)$$
其中,(gcd(a_n,p_n)=1)。
令(k=prod_{i=1}^{n}k_i),则$$k*a^{x_n}equiv b_n(mod p_n)$$
直接使用(BSGS)算法求解该方程,则原方程的解(x=x_n+n)。
(Code:)
inline long long ExBSGS(long long a,long long b,long long p)
{
if (b==1)return 0;
long long cnt=0,d,k=1;
while ( ( d=Euclid(a,p) ) ^ 1 )
{
if (b%d)return -1;
b/=d,p/=d,++cnt;
k = k * (a/d) % p;
if (k==b)return cnt;
}
unordered_map < long long , long long > Hash;Hash.clear();
long long t=(long long)sqrt(p)+1,mul=1;
for (int j=0;j<t;j++)
{
Hash[ mul * b % p ] = j;
mul = mul * a % p;
}
for (int i=0;i<=t;i++)
{
long long j = ( Hash.find(k) != Hash.end() ? Hash[k] : -1 );
if (j>=0&&i*t-j+cnt>=0)return i*t-j+cnt;
k = k * mul % p;
}
return -1;
}
<后记>