NO IMAGE

Note:由於網上很多說一套寫一套的資料,讓我痛苦了整整兩天。在結合了一篇優秀的參考資料後,終於研究得差不多了。我覺得不把這篇總結放出來都有點對不起大家,所以特地把這段總結放上來。其實是希望大家能批評指正一下。

大步小步演算法

1. 問題

​ 求解關於 xx 的方程 ax≡b(modp)a^x \equiv b \; (mod \; p),其中 aa 與 pp 互質。可能無解。

2. 暴力方法

​ O(ϕ(n))O(\phi(n)) 列舉。

3. 大步小步演算法(BSGS)
方法

​ 我們可以設 x=At Bx = At B,tt 是我們自己設的一個引數。其中 0≤B<t0 \le B < t。這樣就變成求 AA 和 BB 了。

aAt B≡b(modp)

a^{At B} \equiv b \; (mod \; p)

​ 由於 aa 與 pp 互質,因此可以方便地求逆元,容易得到:

aAt≡b×a−B(modp)

a^{At} \equiv b × a^{-B} \; (mod \; p)

​ 因此,對於所有的 b×a−Bb × a^{-B},我們可以存進一個雜湊表。然後列舉左邊的 AA,算出 aAta_{At},看看是否在雜湊表中存在。

​ 很明顯引數 tt 設成 ⌈ϕ(c)−−−−√⌉\left \lceil \sqrt {\phi(c)} \right \rceil 是最優的。

​ 這樣預處理的時間複雜度為 O(n√logn)O(\sqrt n \log n),單次查詢的時間複雜度為 O(n√logn)O(\sqrt n \log n)。

4. 擴充套件大步小步演算法(BSGSEx)
①問題

​ 求解關於 xx 的方程 ax≡b(modc)a^x \equiv b \; (mod \; c),其中 aa 與 cc 互質。

②方法

​ 考慮把方程轉換成 gcd(a,c)=1\gcd(a, c) = 1 的形式。

​ 轉換原方程後可以得到這樣一個等價的方程(想想擴充套件歐幾里得):

ax yc=b,y∈Z

a^x yc = b, y \in \mathbb{Z}

​ 由裴蜀定理,若 (g=gcd(a,m))∤b(g = \gcd(a, m)) \nmid b,方程一定無解。

​ 當 g∣bg \mid b 時,在方程左右兩邊同時除以 gg,得到:

agax−1 cgy=bg

\frac {a} {g} a^{x – 1} \frac {c} {g} y = \frac {b} {g}

​ 相當於得到了模方程=

agax−1≡bg(modcg)

\frac {a}{g} a^{x – 1} \equiv \frac {b}{g} \; (mod \; \frac {c}{g})

注意,這裡 gg 不能乘過去了。

​ 令 c′=cgc’ = \frac {c} {g},b′=bg(ag)−1b’ = \frac {b}{g}(\frac {a}{g})^{-1},得新的方程:

ax′≡b′(modc′)

a^{x’} \equiv b’ \; (mod \; c’)

​ 可知 x=x′ 1x = x’ 1。

​ 由於 aa 和 cc 只有 cc 變了,所以可能會有多次過程,不斷遞迴,直到可以用大步小步演算法即可。特別地,當 b′b’ 在某一時刻為 1 時,我們實際上已經得到了唯一解 x′=0x’ = 0。

③實際操作

​ 我們可以用一些奇技淫巧避免求逆元。

(1)對於 BSGS

​ 我們可以設 x=At−Bx = A t – B。其中 0≤B<t0 \le B < t。這樣就變成了:

aAt−B≡b(modp)aAt≡b×aB(modp)

a^{At – B} \equiv b \; (mod \; p)
\\
a^{At} \equiv b × a^{B} \; (mod \; p)

​ 可以避免求逆元。但是從上式推導到下式是要用到逆元的,因此逆元必須存在。

​ tt 要從 1 列舉到 threshold 1threshold 1。

(2)對於 BSGSEx

agax−1≡bg(modcg)

\frac {a}{g} a^{x – 1} \equiv \frac {b}{g} \; (mod \; \frac {c}{g})

​ 由於上式的存在,若要遞迴進行 BSGSEx,我們難以避免求逆元。所以我們可以將整個過程改成非遞迴的。這就要求我們記錄一些東西。

​ 設原式為 ax≡b(modc)a^x \equiv b \; (mod \; c)。每進入一層新的遞迴,我們都會改變 xx,bb 和 cc。對於 xx,我們只會進行減一的操作,所以這個很好維護。對於 cc,我們也只需要讓 cc 不斷除以 (a,c)(a, c) 即可,因為我們只會最後用 cc。稍微複雜點的是 bb,bb 除了要除以 (a,c)(a, c),還要乘以 ag\frac {a} {g} 的逆元。由於每次的操作都是類似的,我們可以儲存每次的 ag\frac {a} {g} 的積(而不是逆元的積)。

​ 因此我們最後的公式看上去是這樣的:(兩個 recrec 代表額外記錄的值)

ax−rec1≡(rec2)−1b′(modc′)

a^{x – rec_1} \equiv (rec_2)^{-1} b’ \; (mod \; c’)

​ rec2rec_2 始終是模 c′c’ 意義下的。不難發現,當 rec2=b‘rec_2 = b‘ 時,x−rec1x – rec_1 有值 00。

​ 可以變形得:

(rec2)ax−rec1≡b′(modc′)

(rec_2)a^{x – rec_1} \equiv b’ \; (mod \; c’)

​ 然後就可以套用大步小步演算法了。設 x′=x−rec1=At−Bx’ = x – rec_1 = At – B。

(rec2)aAt−B≡b′(modc′)(rec2)aAt≡b′×aB(modc′)

(rec_2)a^{At – B} \equiv b’ \; (mod \; c’)
\\
(rec_2)a^{At} \equiv b’ × a^B \; (mod \; c’)

​ 列舉 AA,計算 (rec2)aAt(rec_2)a^{At},在儲存了 b′×aBb’ × a^B 的雜湊表中查詢即可。

5. 參考程式碼
//c   11 is needed
INT BSGSEx(INT a, INT b, INT c) //-1 for no solution
{
a %= c;
b %= c;
if (b == 1)
return 0;
INT count_ = 0; //rec_1
INT base = 1; //rec_2
for (INT g = gcd(a, c); g != 1; g = gcd(a, c))
{
if (b % g) return -1; //exit_1
b /= g;
c /= g;
base = base * (a / g) % c;
count_  ;
if (b == base) return count_; //exit_2
}
INT threshold = std::ceil(std::sqrt(c));
std::unordered_map<INT, INT> table;
INT mul = 1;
for (int i = 0; i < threshold; i  , mul = mul * a % c)
table[mul * b % c] = i;
INT a_t = base * mul % c; //base = rec_2, mul = a ^ threshold
for (int i = 1; i <= threshold   1; i  , a_t = a_t * mul % c) //note the range of i
{
if (table.count(a_t))
return i * threshold - table[a_t]   count_;
}
return -1;
}

裸題傳送門(SPOJ)