形如 az≡b mod pa^z\equiv b \ mod \ paz≡b mod p a,b,pa,b,pa,b,p 已知求 zzz.
如果 gcd(a,p)=1gcd(a,p)=1gcd(a,p)=1
我们令z=xp+y  (x,y<p)z=x \sqrt p +y\; (x,y<\sqrt p)z=xp+y(x,y<p)
az≡b⇒axp+y≡b⇒axp≡b⋅(a−1)ya^z\equiv b \quad \Rightarrow a^{x\sqrt p + y}\equiv b \quad \Rightarrow a^{x\sqrt p}\equiv b\cdot (a^{-1})^yaz≡b⇒axp+y≡b⇒axp≡b⋅(a−1)y
∵ gcd(a,p)=1\because \ gcd(a,p)=1∵ gcd(a,p)=1 我们可以使用扩展欧几里得得到 a−1a^{-1}a−1 .
对于右边 y<py<\sqrt py<p 我们预处理 (a−1)i (0≤i<p)(a^{-1})^i\ (0\le i< \sqrt p)(a−1)i (0≤i<p) 将结果存放在一个表中.
我们枚举左边的 x (0≤x<p)x \ (0\leq x < \sqrt p)x (0≤x<p) 在表中查找即可.
处理 (b⋅a−1)i,axp(b\cdot a^{-1})^i,a^{x\sqrt p}(b⋅a−1)i,axp 我们可以递推得到,复杂度O(p)O(\sqrt p)O(p) . 注意 q\sqrt qq 向上取整
关键在于存储的方式以及查找的复杂度.bitset ?? .unoder_map??. hash?? .
如果 gcd(a,p)≠1gcd(a,p)\neq 1 gcd(a,p)̸=1
az=a⋅az−1≡b⇒a⋅az−1+p⋅y=ba^z=a\cdot a^{z-1}\equiv b\Rightarrow a\cdot a^{z-1} + p\cdot y=baz=a⋅az−1≡b⇒a⋅az−1+p⋅y=b
令 X=az−1X = a^{z-1}X=az−1
⇒a⋅X+p⋅y=b\Rightarrow a\cdot X+p\cdot y=b⇒a⋅X+p⋅y=b 等式两边同时除以 gcd(a,p)gcd(a,p)gcd(a,p) 如果 b%gcd(a,p)≠1b\%gcd(a,p)\neq 1b%gcd(a,p)̸=1无解.
(这里我们假定 XXX 可以取遍 [0,p−1][0,p-1][0,p−1] ).
那么我们得到了一个新式子 a′⋅X+p′⋅y=b′ (a′=a/gcd(a,p),p′=p/gcd(a,p),b′=b/gcd(a,p)a'\cdot X+p'\cdot y=b' \ ( a'=a/gcd(a,p),p'=p/gcd(a,p),b'=b/gcd(a,p)a′⋅X+p′⋅y=b′ (a′=a/gcd(a,p),p′=p/gcd(a,p),b′=b/gcd(a,p)
我们 a′⋅X+p′⋅y=b′⇒a′⋅a⋅az−2≡b′ mod p′a'\cdot X+p'\cdot y=b' \Rightarrow a'\cdot a \cdot a^{z-2} \equiv b' \ mod \ p'a′⋅X+p′⋅y=b′⇒a′⋅a⋅az−2≡b′ mod p′ 令 X=a′⋅az−2X = a'\cdot a^{z-2}X=a′⋅az−2
⇒a⋅X+p′⋅y=b′\Rightarrow a\cdot X+p'\cdot y=b'⇒a⋅X+p′⋅y=b′ 继续迭代到无解或 gcd(a,p′)=1gcd(a,p')=1gcd(a,p′)=1
所以最后的式子可以写成 A⋅az−x=B mod PA \cdot a^{z-x} = B \ mod \ P A⋅az−x=B mod P , gcd(A,P)=1gcd(A,P)=1gcd(A,P)=1
az−x=B⋅A−1 mod Pa^{z-x} = B\cdot A^{-1} \ mod \ Paz−x=B⋅A−1 mod P 利用 gcd(a,P)=1gcd(a,P)=1gcd(a,P)=1 时的情况可以解得 z−xz-xz−x 答案加上 xxx 即可.
这里还存在一个问题,当答案 z<xz<xz<x 时需要特判,由于迭代次数只有几十次 也就是说 xxx 不会很大.暴力检查即可.
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
const int MAXN = 5;
using ll = long long;
ll powmod(ll x,ll n,ll mod) {
ll ret = 1;
for(;n>0;n>>=1,x=x*x%mod) if(n&1) ret = ret*x%mod;
return ret;
}
ll modI(ll a, int m) {
return a<=1?1:(1-modI(m%a, a)*m)/a+m;
}
int BSGS(ll a,ll b,ll p) // a^x=b % p;
{
a %= p,b %= p;
ll a1=1,cnt=0,g;
for(int i=0,res=1;i<100;++i) {
if(res==b) return i;
res = res*a%p;
}
while((g = __gcd(a,p))!=1){
if(b%g) return -1;
b/=g,p/=g;
a1 = a1*(a/g)%p;cnt++;
}
b = modI(a1,p)*b%p;
unordered_map<int,int> mp;
int block = sqrt(p);if(block*block!=p) block++;
ll inv = modI(a,p),res = b;
for(int i=0;i<block;++i) {
if(mp.count(res)==0) mp[res]=i;
res = res*inv%p;
}
res = 1,inv = powmod(a,block,p);
for(int i=0;i<block;++i) {
if(mp.count(res)) return i*block+mp[res]+cnt;
res = res*inv%p;
}
return -1;
}