Menu

数论

发表于 01 Feb 2019 总计 12219字 阅读时间 41分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

辗转相除法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

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$。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

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类型是int的时候小心 a*b 爆精度
	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;
	}
	*/
};

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

使用示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

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;
	}
};

二次剩余

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

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中。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

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];
			}
		}
	}
};

直接求欧拉函数

1
2
3
4
5
6
7
8
9
10
11

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$）：

\[\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)\\ =\sum_d d\sum_{i=1}^n\sum_{j=1}^m[\gcd (i,j)=d]\\ =\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}[\gcd (i,j)=1]\\ =\sum_{d}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac md\rfloor}\sum_{k\mid i,k\mid j}\mu(k)\\ =\sum_d d\sum_k\mu(k)\sum_{k\mid i}^{\lfloor\frac nd\rfloor}\sum_{k\mid j}^{\lfloor\frac md\rfloor}\\ =\sum_{d}d\sum_k\mu(k)\lfloor\frac n{kd}\rfloor\lfloor\frac m{kd}\rfloor\\ =\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\sum_{k\mid T}\frac Tk\mu(k)\\ =\sum_{T}\lfloor\frac nT\rfloor\lfloor\frac mT\rfloor\varphi(T)\]

$\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，看情况选用快速乘法。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

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 时交换的对应位置（即保存的是置换）。

1
2
3
4
5
6
7
8
9

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

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$。

使用示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

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 + 1 r k g，mod(r * 2 ^ k + 1)的原根
3 1 1 2
5 1 2 2
17 1 4 3
97 3 5 5
193 3 6 5
257 1 8 3
7681 15 9 17
12289 3 12 11
40961 5 13 3
65537 1 16 3
786433 3 18 10
5767169 11 19 3
7340033 7 20 3
23068673 11 21 3
104857601 25 22 3
167772161 5 25 3
469762049 7 26 3
998244353 119 23 3
1004535809 479 21 3
2013265921 15 27 31
2281701377 17 27 3
3221225473 3 30 5
75161927681 35 31 3
77309411329 9 33 7
206158430209 3 36 22
2061584302081 15 37 7
2748779069441 5 39 3
6597069766657 3 41 5
39582418599937 9 42 5
79164837199873 9 43 5
263882790666241 15 44 7
1231453023109121 35 45 3
1337006139375617 19 46 3
3799912185593857 27 47 5
4222124650659841 15 48 19
7881299347898369 7 50 6
31525197391593473 7 52 3
180143985094819841 5 55 6
1945555039024054273 27 56 5
4179340454199820289 29 57 3

快速沃尔什变换

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

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

1
2

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

OR

1
2

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

XOR

1
2
3
4
5
6

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$后的结果。