数论 wu-kan

01 Feb 2019 10938字 37分 次 ACM and 算法竞赛模板
如果这篇博客帮助到你，可以请我喝一杯咖啡~
CC BY 4.0 （除特别声明或转载文章外）

辗转相除法

ll gcd(ll a, ll b) { return b ? gcd(b, a % b) : a; } //a、b的最大公约数
ll lcm(ll a, ll b) { return a / gcd(a, b) * b; }	 //a、b的最小公倍数
ll gcd(ll a, ll b, ll &x, ll &y)					 //扩展欧几里得，引用返回a*x+b*y=gcd(a,b)绝对值之和最小的解
{
	if (!a)
		return x = 0, y = 1, b;
	ll d = gcd(b % a, a, y, x);
	return x -= b / a * y, d;
}
ll gcd(ll a, ll b) //无取模求gcd
{
	for (ll t = 1, c, d;;)
	{
		if (a == b)
			return t * a;
		if (a < b)
			swap(a, b);
		if (a & 1)
			c = 0;
		else
			a >>= 1, c = 1;
		if (b & 1)
			d = 0;
		else
			b >>= 1, d = 1;
		if (c && d)
			t <<= 1;
		else if (!c && !d)
			a -= b;
	}
}

裴蜀定理

对任何$a,b\in Z$和它们的最大公约数$d$，关于未知数$x,y$的线性不定方程（称为裴蜀等式）：$ax+by=c$当仅当$d\vert c$，可知有无穷多解。特别地，$ax+by=d$一定有解。

推论

$a,b$互质的充要条件是$ax+by=1$有整数解。

同余系运算

三种乘法的比较，结论是如果可以用__int128就用，否则就用第一种：

求乘法逆元的另外一种方法是用欧拉定理$x^{\phi(m)}\equiv1\pmod m$，x 的逆是$x^{\phi(m)-1}$。特别地，m 为素数时$\phi(m)=m-1$，此时 x 的逆就是pow(x,m-2,m)。

log 函数：m 为素数时求解模方程$a^x\equiv b\pmod m$。设 P 为质数，G 为 P 的原根，则$x^y\equiv b\pmod P$等价于$y\ ind\ x\equiv b\pmod{P-1}$，其中$G\ ind\ x\equiv x\pmod P$。

struct Mod
{
	const ll M, SM;
	Mod(ll M) : M(M), SM(sqrt(M) + 0.5) {}
	ll qadd(ll &a, ll b) const { return a += b, a >= M ? a -= M : a; } //假如a+b<2*M，就不必取模了，取模运算耗时很高
	ll add(ll a, ll b) const { return qadd(a = (a + b) % M, M); }	  //考虑a和b不在同余系内甚至为负数的情况
	ll mul(ll a, ll b) const { return add(a * b, M); }
	ll inv(ll a) const { return pow(a, M - 2); } //要求M为素数，否则return pow(a, phi(M) - 1);
	ll pow(ll a, ll b) const
	{
		ll r = 1;
		/*
		if (b < 0)
			b = -b, a = inv(a);
		*/
		for (a = add(a, M); b; b >>= 1, a = mul(a, a))
			if (b & 1)
				r = mul(r, a);
		return r;
	}
	/*
	ll mul(ll a, ll b) const { return add(a * b, -M * (ll)((long double)a / M * b)); } //long double 有效位数16~18位，模数过大时慎用！
	ll mul(ll a, ll b) const //Head算法，无循环快速计算同余乘法，根据a*b是否爆ll替换a*b%M，需要a<M且b<M，可以调用时手动取模
	{
		ll c = a / SM, d = b / SM;
		a %= SM, b %= SM;
		ll e = add(add(a * d, b * c), c * d / SM * (SM * SM - M));
		return add(add(a * b, e % SM * SM), add(c * d % SM, e / SM) * (SM * SM - M));
	}
	ll mul(ll a, ll b) const //龟速乘
	{
		ll r = 0;
		for (a %= M; b; b >>= 1, qadd(a, a))
			if (b & 1)
				qadd(r, a);
		return r;
	}
	ll inv(ll a) const //模m下a的乘法逆元，不存在返回-1（m为素数时a不为0必有逆元）
	{
		ll x, y, d = gcd(a, M, x, y);
		return d == 1 ? add(x, M) : -1;
	}
	vector<ll> sol(ll a, ll b) const //解同余方程，返回ax=b(mod M)循环节内所有解
	{
		vector<ll> ans;
		ll x, y, d = gcd(a, M, x, y);
		if (b % d)
			return ans;
		ans.push_back(mul((b / d) % (M / d), x));
		for (ll i = 1; i < d; ++i)
			ans.push_back(add(ans.back(), M / d));
		return ans;
	}
	ll log(ll a, ll b) const
	{
		unordered_map<ll, ll> x;
		for (ll i = 0, e = 1; i <= SM; ++i, e = mul(e, a))
			if (!x.count(e))
				x[e] = i;
		for (ll i = 0, v = inv(pow(a, SM)); i <= SM; ++i, b = mul(b, v))
			if (x.count(b))
				return i * SM + x[b];
		return -1;
	}
	*/
};

中国剩余定理解同余方程组

使用示例

ll crt(const vector<pair<ll, ll>> &v) //同余方程组，x%v[i].first==v[i].second，不存在返回-1
{
	ll m = v[0].first, r = v[0].second, c, d, x, y, z;
	for (int i = 1; i < v.size(); ++i)
	{
		if (c = v[i].second - r, d = gcd(m, v[i].first, x, y), c % d)
			return -1;
		gcd(m / d, z = v[i].first / d, x, y), r += c / d * x % z * m, r %= m *= z;
	}
	return r < 0 ? r + m : r;
}
struct Excrt
{
	vector<ll> b2, a2;
	int n;
	void init()
	{
		scanf("%d", &n);
		ll q, qq;
		for (int i = 0; i < n; i++)
			scanf("%lld%lld", &q, &qq), b2.push_back(q), a2.push_back(qq);
	}
	ll get()
	{
		ll M, ans, x = 0, y = 0;
		M = 1, ans = 0;
		for (int i = 0; i < n; i++)
		{
			ll a = M, b = b2[i], c = (a2[i] - ans % b + b) % b; //ax=c(mod b)
			ll gcd = exgcd(a, b, x, y), bg = b / gcd;
			if (c % gcd != 0)
			{
				return -1;
			} //返回-1表示无解
			x = mul(x, c / gcd, bg);
			y = M / gcd * b;
			ans = (x * M + ans) % y; //???k???????
			if (ans < 0)
				ans += y;
			M = y; //M??k?m?lcm
		}
		return ans;
	}
};

二次剩余

struct QuadraticRes
{
	struct numb
	{
		ll x, y;
	};
	numb mul(const numb &a, const numb &b, ll p)
	{
		return (numb){(a.x * b.x % p + a.y * b.y % p * w % p) % p, (a.x * b.y + a.y * b.x) % p};
	}
	ll pow2(numb a, ll b, ll p)
	{
		numb ret = {1, 0}, tem = a;
		while (b)
		{
			if (b & 1)
				ret = mul(ret, tem, p);
			tem = mul(tem, tem, p), b >>= 1;
		}
		return ret.x % p;
	}
	ll get(ll n, ll p)
	{
		if (Mod(p).pow(n, p / 2) == p - 1)
		{
			return -1;
		} //此时无解
		while (1)
		{
			ll a = rand();
			w = (a * a - n + p) % p;
			if (Mod(p).pow(w, p / 2) == p - 1)
			{
				return pow2((numb){a, 1}, p / 2 + 1, p); //还有一个解是p减去这个值
			}
		}
	}
};

欧拉筛

欧拉函数$\phi(n)$是小于 n 的正整数中与 n 互素的数的数目。特别地，规定$\phi(1)=1$，易知 n>2 时都为偶数。

欧拉函数是积性函数，即对任意素数$p,q$满足下列关系：$\phi(pq)=\phi(p)\phi(q)=(p-1)(q-1)$对任何两个互质的正整数$x, m(m\geq2)$有欧拉定理：$x^{\phi(m)}\equiv1\pmod m$当 m 为素数 p 时，此式变为费马小定理：$x^{p-1}\equiv1\pmod p$利用欧拉函数和它本身不同质因数的关系，用筛法$O(N)$预处理某个范围内所有数的欧拉函数值，并求出素数表。同时，利用计算欧拉函数过程中求出的最小素因子 m，可以实现$O(log N)$的素因数分解。

同时求莫比乌斯函数$\mu(n)$，存在mu中。

struct EulerSieve
{
	vector<int> p, m, phi, mu; //素数序列，最小素因子，欧拉函数，莫比乌斯函数
	EulerSieve(int N) : m(N, 0), phi(N, 0), mu(N, 0)
	{
		phi[1] = mu[1] = 1;					 //m[1]=0,m[i]==i可判断i是素数
		for (long long i = 2, k; i < N; ++i) //防i*p[j]爆int
		{
			if (!m[i])
				p.push_back(m[i] = i), phi[i] = i - 1, mu[i] = -1; //i是素数
			for (int j = 0; j < p.size() && (k = i * p[j]) < N; ++j)
			{
				phi[k] = phi[i] * p[j];
				if ((m[k] = p[j]) == m[i])
				{
					mu[k] = 0;
					break;
				}
				phi[k] -= phi[i];
				mu[k] = -mu[i];
			}
		}
	}
};

直接求欧拉函数

ll phi(ll n)
{
	ll phi = n;
	for (ll i = 2; i * i <= n; ++i)
		if (!(n % i))
			for (phi = phi / i * (i - 1); !(n % i);)
				n /= i;
	if (n > 1)
		phi = phi / n * (n - 1);
	return phi;
}

常见数论函数变换

$\sum_{d\vert n}\mu(d)=[n=1]$

$\phi(n)=\sum_{i=1}^n[\gcd(i,n)=1]=\sum_{i=1}^n\sum_{k\mid i,k\mid n}\mu(k)=\sum_{k\mid n}\frac nk\mu(k)$

前缀和

欧拉函数前缀和$S_\phi(n)=\frac{(n+1)n}2-\sum_{d=1}^nS_\phi(\frac{n}{d})$

莫比乌斯函数前缀和$S_\mu(n)=1-\sum_{d=1}^nS_\mu(\frac{n}{d})$

莫比乌斯反演

若$f(n)=\sum_{d\vert n}g(d)$，则$g(n)=\sum_{d\vert n}\mu(d)f(\frac{n}{d})$

若$f(n)=\sum_{i=1}^nt(i)g(\frac{n}{i})$，则$g(n)=\sum_{i=1}^n\mu(i)t(i)f(\frac{n}{i})$（此时代$t(i)=[gcd(n,i)>1]$可得上式）

举例（其中$T=kd$）：

$\varphi(T)$可以使用线性筛预处理处理，我们就可以枚举$T$求上式了，时间复杂度$O(n)$。多组数据$n,m$询问上式，时间复杂度就变成了$O(Tn)$。事实上，$\lfloor\frac{n}{T}\rfloor$是不会轻易变化的，是过了连续的一段后才发生变化的，那么我们就可以计算出这一段的结束位置，对$\varphi$函数作前缀和，就可以直接分块了，这样的时间复杂度是$O(T\sqrt{n})$的。

PollardRho 大数素因子分解

时间复杂度$O(N^{1/4})$，数据多的时候可考虑欧拉筛优化。

注意这里模乘法很容易爆 long long，看情况选用快速乘法。

struct PollardRho
{
	bool isPrime(ll n, int S = 12) //MillerRabin素数测试，S为测试次数，用前S个素数测一遍，S=12可保证unsigned long long范围内无错；n<2请特判
	{
		static ll d, u, t, p[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
		for (d = n - 1; !(d & 1);)
			d >>= 1; //未对0，1做特判！
		Mod mo(n);
		for (int i = 0; i < S; ++i)
		{
			if (!(n % p[i]))
				return n == p[i];
			for (t = mo.pow(p[i], u = d); t != n - 1 && t != 1 && u != n - 1;)
				t = mo.mul(t, t), u <<= 1;
			if (t != n - 1 && !(u & 1))
				return 0;
		}
		return 1;
	}
	void fac(ll n, vector<ll> &factor)
	{
		/*
		if (n < e.m.size()) //欧拉筛预处理优化，可以防止很多小素数因子的情况
		{
			for (; n > 1; n /= e.m[n])
				factor.push_back(e.m[n]);
			return;
		}
		*/
		if (isPrime(n))
			return factor.push_back(n);
		Mod mo(n);
		for (ll c = 1;; ++c)
			for (ll i = 0, k = 1, x = rand() % (n - 1) + 1, y, p;;)
			{
				if (++i == k)
					y = x, k <<= 1;
				if (x = mo.add(mo.mul(x, x), c), p = __gcd(abs(x - y), n), p == n)
					break;
				if (p > 1)
					return fac(p, factor), fac(n / p, factor);
			}
	}
};

快速变换

蝴蝶变换（雷德变换）

保存 FTT 和 FNTT 时交换的对应位置（即保存的是置换）。

struct Rader : vector<int>
{
	Rader(int n) : vector<int>(1 << int(ceil(log2(n))))
	{
		for (int i = at(0) = 0; i < size(); ++i)
			if (at(i) = at(i >> 1) >> 1, i & 1)
				at(i) += size() >> 1;
	}
};

快速傅里叶变换

使用示例，高精度乘法。

任意模数 FFT

struct FFT : Rader
{
	vector<complex<lf>> w;
	FFT(int n) : Rader(n), w(size())
	{
		for (int i = 0; i < size(); ++i)
			w[i] = polar(1.0, 2 * M_PI * i / size());
	}
	void fft(vector<complex<lf>> &x) const
	{
		for (int i = 0; i < x.size(); ++i)
			if (i < at(i))
				std::swap(x[i], x[at(i)]);
		for (int i = 1; i < size(); i <<= 1)
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					complex<lf> t = w[size() / (i << 1) * j] * x[k + i];
					x[k + i] = x[k] - t, x[k] += t;
				}
	}
	vector<ll> ask(const vector<ll> &a, const vector<ll> &b) const
	{
		vector<complex<lf>> xa(a.begin(), a.end()), xb(b.begin(), b.end());
		xa.resize(size()), xb.resize(size()), fft(xa), fft(xb);
		for (int i = 0; i < size(); ++i)
			xa[i] *= xb[i];
		fft(xa);
		vector<ll> ans(size());
		for (int i = 0; i < size(); ++i)
			ans[i] = xa[(size() - i) & (size() - 1)].real() / size() + 0.5;
		return ans;
	}
	vector<ll> askMod(const vector<ll> &a, const vector<ll> &b, ll M) const //任意模数
	{
		vector<complex<lf>> e(size()), c(size());
		for (int i = 0; i < a.size(); ++i)
			e[i].real(a[i] & 0x7fff), e[i].imag(a[i] >> 15);
		for (int i = 0; i < b.size(); ++i)
			c[i].real(b[i] & 0x7fff), c[i].imag(b[i] >> 15);
		fft(e), fft(c);
		vector<complex<lf>> d(c);
		for (int i = 0; i < size(); ++i)
		{
			int fr = (size() - i) & (size() - 1);
			c[i] *= 0.5 * complex<lf>(e[i].X + e[fr].X, e[i].Y - e[fr].Y);
			d[i] *= 0.5 * complex<lf>(e[i].Y + e[fr].Y, e[fr].X - e[i].X);
		}
		fft(c), fft(d);
		vector<ll> ans(size());
		for (int i = 0; i < size(); ++i)
		{
			ll fr = (size() - i) & (size() - 1),
			   p = c[fr].X / size() + 0.5,
			   o = c[fr].Y / size() + 0.5,
			   x = d[fr].X / size() + 0.5,
			   u = d[fr].Y / size() + 0.5;
			ans[i] = (p % M + ((o + x) % M << 15) + (u % M << 30)) % M;
		}
		return ans;
	}
};

快速数论变换

原理和 FFT 相同，解决特殊情况下 FFT 的浮点误差，并且可以在同余系进行变换。

对于形如$m=2^nc+1$的费马素数，记其原根为 g，则旋转因子为$g^{(m-1)/n}$，满足$g^{m-1}=1$且$2^n\mid m-1$。

使用示例

struct FNTT : Rader, Mod
{
	vector<ll> w;
	FNTT(int N, ll M, ll G) : Rader(N), Mod(M), w(size(), pow(G, (M - 1) / size()))
	{
		for (int i = w[0] = 1; i < size(); ++i)
			w[i] = mul(w[i], w[i - 1]);
	}
	void fntt(vector<ll> &x) const
	{
		for (int i = 0; i < size(); ++i)
			if (i < at(i))
				std::swap(x[i], x[at(i)]);
		for (int i = 1, j; i < size(); i <<= 1)
			for (int j = 0; j < i; ++j)
				for (int k = j; k < size(); k += i << 1)
				{
					ll t = mul(w[size() / (i << 1) * j], x[k + i]);
					qadd(x[k + i] = x[k], M - t), qadd(x[k], t);
				}
	}
	vector<ll> ask(vector<ll> a, vector<ll> b) const
	{
		a.resize(size()), b.resize(size()), fntt(a), fntt(b);
		for (int i = 0; i < size(); ++i)
			a[i] = mul(a[i], b[i]);
		fntt(a), reverse(a.begin() + 1, a.end());
		ll u = inv(size());
		for (int i = 0; i < size(); ++i)
			a[i] = mul(a[i], u);
		return a;
	}
};

常见素数的原根

r * 2 ^ k + 1rkg，mod(r * 2 ^ k + 1)的原根
3112
5122
17143
97355
193365
257183
768115917
1228931211
409615133
655371163
78643331810
576716911193
73400337203
2306867311213
10485760125223
1677721615253
4697620497263
998244353119233
1004535809479213
2013265921152731
228170137717273
32212254733305
7516192768135313
773094113299337
20615843020933622
206158430208115377
27487790694415393
65970697666573415
395824185999379425
791648371998739435
26388279066624115447
123145302310912135453
133700613937561719463
379991218559385727475
4222124650659841154819
78812993478983697506
315251973915934737523
1801439850948198415556
194555503902405427327565
417934045419982028929573

快速沃尔什变换

如果要在同余系中进行运算，则下面代码需要修改。

void fwt(vector<ll> &x, void f(ll &l, ll &r))
{
	for (int i = 1; i < x.size(); i <<= 1)
		for (int j = 0; j < i; ++j)
			for (int k = j; k < x.size(); k += i << 1)
				f(x[k], x[k + i]);
}
void fwt(ll *b, ll *e, void f(ll &l, ll &r)) //再给一个递归二分的代码便于理解
{
	if (e - b < 2)
		return;
	ll *m = b + (e - b) / 2;
	fwt(b, m, f), fwt(m, e, f);
	while (m < e)
		f(*(b++), *(m++));
}

AND

void tf(ll &l, ll &r) { l += r; }
void utf(ll &l, ll &r) { l -= r; }

OR

void tf(ll &l, ll &r) { r += l; }
void utf(ll &l, ll &r) { r -= l; }

XOR

void tf(ll &l, ll &r)
{
	ll tl = l + r, tr = l - r;
	l = tl, r = tr;
}
void utf(ll &l, ll &r) { tf(l, r), l >>= 1, r >>= 1; }

XNOR、NAND、NOR

直接用异或运算、与运算、或运算的方法求出来，然后将互反的两位交换即可。

Pell 方程

形如$x^2-Dy^2=1$（D 为任意正整数）的方程称为佩尔方程，必有最小正整数解$(x_0,y_0)$，用

$x_n=x_0x_{n-1}+Dy_0y_{n-1},y_n=y_0x_{n-1}+x_0y_{n-1}$

可递推方程的第 n 小整数解（可用矩阵快速幂求），同时还有

$2x_0x_n=x_{n-1}+x_{n+1},2x_0y_n=y_{n-1}+y_{n+1}$

Bertrand 猜想

$\forall n>3,\exist n<p<n\times 2$其中 n 为整数，p 为质数。

威尔逊定理

当且仅当 p 为素数时：$(p-1)!\equiv -1\pmod p$，

Jacobi’s Four Square Theorem

设$a^2+b^2+c^2+d^2=n$的自然数解个数为$r4(n)$，$d(n)$为$n$的约数和，由 Jacobi’s Four Square Theorem 可知，若$n$是奇数，则$r4(n)=8d(n)$，否则$r4(n)=24d(k)$，$k$是$n$去除所有$2$后的结果。