数论 · dfydada

辗转相除法

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都已经在同余系内，就不必取模了，取模运算耗时很高
	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;
		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;
}

欧拉筛

欧拉函数$\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$）：

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

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

快速傅里叶变换

使用示例，高精度乘法。

struct FFT : Rader
{
	vector<complex<lf>> w;
	FFT(int n) : Rader(n), w(size(), polar(1.0, 2 * M_PI / size()))
	{
		w[0] = 1;
		for (int i = 1; i < size(); ++i)
			w[i] *= w[i - 1];
	}
	vector<complex<lf>> fft(const vector<complex<lf>> &a) const
	{
		vector<complex<lf>> x(size());
		for (int i = 0; i < a.size(); ++i)
			x[at(i)] = a[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;
				}
		return x;
	}
	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 = fft(xa), xb = fft(xb);
		for (int i = 0; i < size(); ++i)
			xa[i] *= xb[i];
		vector<ll> ans(size());
		xa = fft(xa), ans[0] = xa[0].real() / size() + 0.5;
		for (int i = 1; i < size(); ++i)
			ans[i] = xa[size() - i].real() / size() + 0.5;
		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
{
	ll G;
	vector<ll> w;
	FNTT(int N, ll M, ll G) : Rader(N), Mod(M), G(G), w(size(), pow(G, (M - 1) / size()))
	{
		for (int i = w[0] = 1; i < size(); ++i)
			w[i] = mul(w[i], w[i - 1]);
	}
	vector<ll> fntt(const vector<ll> &a) const
	{
		vector<ll> x(size());
		for (int i = 0; i < a.size(); ++i)
			x[at(i)] = a[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]);
					x[k + i] = qadd(x[k], M - t), x[k] = qadd(x[k], t);
				}
		return x;
	}
	vector<ll> ask(vector<ll> a, vector<ll> b) const
	{
		a = fntt(a), b = fntt(b);
		for (int i = 0; i < size(); ++i)
			a[i] = mul(a[i], b[i]);
		a = 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;
	}
};

快速沃尔什变换

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

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

数论 dfydada