給定整數 \(a,b,p\) 且 \(a,p\) 互質，請求出高次同餘方程

\[a^x\equiv b\pmod p \]

的非負整數解.

首先，

\[a^0\equiv 1\pmod p \]

又由尤拉定理知

\[a^{\varphi(p)}\equiv 1\pmod p \]

所以出現了週期.

所以 \(x\) 在區間 \([0,\varphi(p))\) 中一定有，否則就是無解.

那麼暴力就是遍歷 \(0\sim \varphi(p)-1\)​ 去試每一個，時間複雜度為 \(\mathcal{O}(\varphi(p))\)​.

---

\(\rm Baby\ Step,Giant\ Step\)，簡稱 \(\rm BSGS\)

演算法，俗稱大小步演算法.

當 \(a,p\) 不互質時就得用 \(\rm EXBSGS\) 了（然而我太菜了不會

把 \(x\) 表示成帶餘除法的樣子：\(x=it+j\)，其中 \(t\) 是某個整數（取值待會討論），那麼有 \(0\le i\le\left\lfloor\dfrac{\varphi(p)}{t} \right\rfloor,0\le j\le t-1\).

\[a^x\equiv b\pmod p\\ a^{it+j}\equiv b\pmod p\\ a^{it}\equiv b\cdot a^{-j}\pmod p\\ (a^t)^i\equiv b\cdot a^{-j}\pmod p \]

\(\rm Baby\ Step\)

：\(\forall j\in [0,t-1]\)，將 \((b\cdot a^{-j})\bmod p\) 扔進 \(\rm Hash\) 表.

\(\rm Giant\ Step\)：\(\forall i\in\left[0,\left\lfloor \dfrac{\varphi(p)}{t}\right\rfloor\right]\)，在 \(\rm Hash\) 表中找是否有 \((a^t)^i\bmod p\)，若有，則 \(x=it+j\) 就是其中一個解.

\(\rm Baby\ Step\) 的時間複雜度為 \(\mathcal{O}(t\log t)\)，\(\rm Giant\ Step\)

的時間複雜度為 \(\mathcal{O}\left(\left\lfloor\dfrac{\varphi(p)}{t} \right\rfloor\log t\right)\)，所以總的時間複雜度就是 \(\mathcal{O}\left(\max\left(t,\left\lfloor\dfrac{\varphi(p)}{t} \right\rfloor\right)\log t\right)\)，這個很像分塊的時間，當取 \(t=\left\lfloor\sqrt{\varphi(p)} \right\rfloor\) 時最優.

但是有一個問題，上式右邊是 \(b\cdot a^{-j}\)，導致需要處理逆元，肥腸麻煩.

所以我們改為令 \(x=it-j\).

---

\(x=it-j\)，那麼有 $ 0\le i\le \left\lceil\dfrac{\varphi(p)}{t}\right\rceil$，\(0\le j\le t-1\).

\[a^x\equiv b\pmod p\\ a^{it-j}\equiv b\pmod p\\ (a^t)^i\equiv b\cdot a^j\pmod p \]

\(\rm Baby\ Step\)：\(\forall j\in [0,t-1]\)，將 \((b\cdot a^{j})\bmod p\) 扔進 \(\rm Hash\) 表.

\(\rm Giant\ Step\)：\(\forall i\in\left[0,\left\lceil \dfrac{\varphi(p)}{t}\right\rceil\right]\)，在 \(\rm Hash\) 表中找是否有 \((a^t)^i\bmod p\)，若有且 \(it-j\ge 0\)，則 \(x=it-j\) 就是其中一個解.

總的時間複雜度是 \(\mathcal{O}\left(\max\left(t,\left\lceil\dfrac{\varphi(p)}{t} \right\rceil\right)\log t\right)\)，當取 \(t=\left\lceil\sqrt{\varphi(p)} \right\rceil\) 時最優.

P3846 [TJOI2007] 可愛的質數/【模板】BSGS

對於本題，要求出最小解，因為 \(0\le j\le t-1\)，所以要使 \(x\) 最小隻需使 \(i\) 最小，按順序遍歷，第一次找到就返回.

又因為 \(p\) 為質數，\(\varphi(p)=p-1\)，則這裡的 \(t\) 都取 \(\left\lceil\sqrt{p}\right\rceil\) 即可.

此時，\(i\) 從 \(0\) 迴圈到 \(t\) 即可.

\(\text{Code}\)

#include <iostream>
#include <cstdio>
#include <map>
#include <cmath>
using namespace std;

int BSGS(int a, int b, int p)
{
	b %= p;
	map<int, int> hash;
	int t = ceil(sqrt(p)), aj = 1;
	for (int j = 0; j < t; j++)
	{
		int val = (long long)b * aj % p;
		hash[val] = j;
		aj = (long long)aj * a % p;
	}
	a = aj; //a^t
	if (!a)
	{
		if (!b)
		{
			return 1;
		}
		else
		{
			return -1;
		}
	}
	int ai = 1;
	for (int i = 0; i <= t; i++)
	{
		int j = hash[ai];
		if (j > 0 && i * t - j >= 0) //這裡必須保證結果非負
		{
			return i * t - j;
		}
		ai = (long long)ai * a % p;
	}
	return -1;
}

int main()
{
	int a, b, p;
	scanf("%d%d%d", &p, &a, &b);
	int res = BSGS(a, b, p);
	if (res == -1)
	{
		puts("no solution");
	}
	else
	{
		printf("%d\n", res);
	}
	return 0;
}

---

矩陣 \(\rm BSGS\)

【BZOJ4128】Matrix

給定矩陣 \(A,B\) 和模數 \(P\)，求最小的 \(x\)，滿足 \(A^x\equiv B\pmod P\).

用矩陣乘法 \(+\) 矩陣 \(\rm Hash\) 即可.

由於常數巨大，需要一些玄學優化：將原來的

for (int i = 1; i <= n; i++)
{
	for (int j = 1; j <= n; j++)
	{
		for (int k = 1; k <= n; k++)
		{
			res.a[i][j] = (res.a[i][j] + a[i][k] * x.a[k][j]) % p;
		}
	}
}

改成

for (int k = 1; k <= n; k++)
{
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			res.a[i][j] = (res.a[i][j] + a[i][k] * x.a[k][j]) % p;
		}
	}
}

這樣就可以剪個枝

for (int k = 1; k <= n; k++)
{
	for (int i = 1; i <= n; i++)
	{
		if (!a[i][k])
		{
			continue;
		}
		for (int j = 1; j <= n; j++)
		{
			res.a[i][j] = (res.a[i][j] + a[i][k] * x.a[k][j]) % p;
		}
	}
}

\(\text{Code}\)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <map>
#include <cmath>
using namespace std;

const int MAXN = 75;
const int BASE = 233;
const int MOD = 19260817;

int n, p;

struct matrix
{
	int a[MAXN][MAXN];
	void init()
	{
		memset(a, 0, sizeof(a));
	}
	void build()
	{
		for (int i = 1; i <= n; i++)
		{
			a[i][i] = 1;
		}
	}
	matrix operator *(const matrix &x)const
	{
		matrix res;
		res.init();
		for (int k = 1; k <= n; k++)
		{
			for (int i = 1; i <= n; i++)
			{
				if (!a[i][k])
				{
					continue;
				}
				for (int j = 1; j <= n; j++)
				{
					res.a[i][j] = (res.a[i][j] + a[i][k] * x.a[k][j]) % p;
				}
			}
		}
		return res;
	}
	int hash()
	{
		int ans = 0;
		for (int i = 1; i <= n; i++)
		{
			for (int j = 1; j <= n; j++)
			{
				ans = (ans * BASE + a[i][j]) % MOD;
			}
		}
		return ans;
	}
};

int BSGS(matrix a, matrix b)
{
	map<int, int> hash;
	int t = ceil(sqrt(p));
	matrix aj;
	aj.init(), aj.build();
	for (int j = 0; j < t; j++)
	{
		matrix val = b * aj;
		hash[val.hash()] = j;
		aj = aj * a;
	}
	a = aj;
	matrix ai;
	ai.init(), ai.build();
	for (int i = 0; i <= t; i++)
	{
		int j = hash[ai.hash()];
		if (j > 0 && i * t - j >= 0)
		{
			return i * t - j;
		}
		ai = ai * a;
	}
}

int main()
{
	scanf("%d%d", &n, &p);
	matrix a;
	a.init();
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			scanf("%d", &a.a[i][j]);
		}
	}
	matrix b;
	b.init();
	for (int i = 1; i <= n; i++)
	{
		for (int j = 1; j <= n; j++)
		{
			scanf("%d", &b.a[i][j]);
		}
	}
	printf("%d\n", BSGS(a, b));
	return 0;
}