[模板]数学

Posted on Edited on In

[TOC]

答案不取模，亲人两行泪

212370440130137957

(↑这是一个很大的质数)

19260817

(↑这是一颗很大的子弹)

判断同余时用\((a - b) \mod m== 0\)

cout不使用科学计数法输出，设置精度: 

1
2
cout.setf(ios::fixed,ios::floatfield);
cout.precision(2);

# 杂项
## 神秘结论
\[ \sum_{i = 1}^n (d(i))^2 \ll n \log^3 n (d为因子个数) \]
## 大整数
(我也不知道从哪copy)的 
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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
struct bign{
    int d[maxn], len;
 
	void clean() { while(len > 1 && !d[len-1]) len--; }
 
    bign() 			{ memset(d, 0, sizeof(d)); len = 1; }
    bign(int num) 	{ *this = num; } 
	bign(char* num) { *this = num; }
    bign operator = (const char* num){
        memset(d, 0, sizeof(d)); len = strlen(num);
        for(int i = 0; i < len; i++) d[i] = num[len-1-i] - '0';
        clean();
		return *this;
    }
    bign operator = (int num){
        char s[20]; sprintf(s, "%d", num);
        *this = s;
		return *this;
    }
 
    bign operator + (const bign& b){
        bign c = *this; int i;
        for (i = 0; i < b.len; i++){
        	c.d[i] += b.d[i];
        	if (c.d[i] > 9) c.d[i]%=10, c.d[i+1]++;
		}
		while (c.d[i] > 9) c.d[i++]%=10, c.d[i]++;
		c.len = max(len, b.len);
		if (c.d[i] && c.len <= i) c.len = i+1;
        return c;
    }
    bign operator - (const bign& b){
        bign c = *this; int i;
        for (i = 0; i < b.len; i++){
        	c.d[i] -= b.d[i];
        	if (c.d[i] < 0) c.d[i]+=10, c.d[i+1]--;
		}
		while (c.d[i] < 0) c.d[i++]+=10, c.d[i]--;
		c.clean();
		return c;
    }
    bign operator * (const bign& b)const{
        int i, j; bign c; c.len = len + b.len; 
        for(j = 0; j < b.len; j++) for(i = 0; i < len; i++) 
			c.d[i+j] += d[i] * b.d[j];
        for(i = 0; i < c.len-1; i++)
            c.d[i+1] += c.d[i]/10, c.d[i] %= 10;
        c.clean();
		return c;
    }
    bign operator / (const bign& b){
    	int i, j;
		bign c = *this, a = 0;
    	for (i = len - 1; i >= 0; i--)
    	{
    		a = a*10 + d[i];
    		for (j = 0; j < 10; j++) if (a < b*(j+1)) break;
    		c.d[i] = j;
    		a = a - b*j;
    	}
    	c.clean();
    	return c;
    }
    bign operator % (const bign& b){
    	int i, j;
		bign a = 0;
    	for (i = len - 1; i >= 0; i--)
    	{
    		a = a*10 + d[i];
    		for (j = 0; j < 10; j++) if (a < b*(j+1)) break;
    		a = a - b*j;
    	}
    	return a;
    }
	bign operator += (const bign& b){
        *this = *this + b;
        return *this;
    }
 
    bool operator <(const bign& b) const{
        if(len != b.len) return len < b.len;
        for(int i = len-1; i >= 0; i--)
            if(d[i] != b.d[i]) return d[i] < b.d[i];
        return false;
    }
    bool operator >(const bign& b) const{return b < *this;}
    bool operator<=(const bign& b) const{return !(b < *this);}
    bool operator>=(const bign& b) const{return !(*this < b);}
    bool operator!=(const bign& b) const{return b < *this || *this < b;}
    bool operator==(const bign& b) const{return !(b < *this) && !(b > *this);}
 
    string str() const{
        char s[maxn]={};
        for(int i = 0; i < len; i++) s[len-1-i] = d[i]+'0';
        return s;
    }
};
 
istream& operator >> (istream& in, bign& x)
{
    string s;
    in >> s;
    x = s.c_str();
    return in;
}
 
ostream& operator << (ostream& out, const bign& x)
{
    out << x.str();
    return out;
}
使用例: 
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
bign GCD(bign a,bign b)
{
	// cout<<"@@@ "<<a<<endl<<"### "<<b<<endl;;
	bign k=0,c;
	while(b>k)
	{
		c=a%b;
		a=b;
		b=c;
	}
	return a;
}
int T;
int vis[110];
vector <bign> prime;
void GetPrime()
{
	for(long long i=2;;i+=1)
	{
		int flag=1;
		for(long long j=2;j<i;j+=1)
			if(i%j==0) flag=0;
		if(flag)
		{
			bign k=i;
			prime.push_back(k);
		}
		if((int)prime.size()>100) break;
	}
}
int main()
{
	cin>>T;
	GetPrime();
	while(T--)
	{
		bign n;
		cin>>n;
		bign a=1,b=1,k=1,sum=1;
		int pos=-1;
		for(int i=0;i<100;i++)
		{
			if(sum*prime[i]<=n) sum=sum*prime[i],pos=i;
			else break;
		}
 
		for(int i=0;i<=pos;i++)
		{
			
			b=b*prime[i];
			a=a*(prime[i]+k);
			
			// cout<<"@@@ "<<s<<endl;
			
		}
		bign s=GCD(a,b);
		b=b/s,a=a/s;
		cout<<b<<"/"<<a<<endl;
	}
	return 0;
 
}
## 随机数发生器
注意不要用Windows下的\(rand()\) 
1
2
3
mt19937 Rand(seed);
mt19937 Rand(chrono::steady_clock::now().time_since_epoch().count());//cf最好用这个防止被hack
Rand();//返回一个随机数
### 均匀分布 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
// 包含random头文件
#include<random>
#include<iostream>
int main() {
	//定义均匀分布对象，均匀分布区间（a,b)为(-10，10)
	std::uniform_int_distribution<int> uid{ -10,10 };
	//均匀分布区间(a,b)
	std::cout << uid.a() <<" "<< uid.b() << std::endl;
	//定义随机数种子
	std::random_device rd;
	//定义默认随机数生成器
	std::default_random_engine dre{ rd() };
	std::cout << "五个随机数如下" << std::endl;
	for (int i = 0; i < 5; ++i)
		std::cout << uid(dre) << std::endl;
	return 0;

}
## 计时 
1
2
3
4
clock_t start = clock();
// do something...
clock_t end   = clock();
cout << "花费了" << (double)(end - start) / CLOCKS_PER_SEC << "秒" << endl;
# 数论
## 同余式
同余属于等价关系, 满足自反、对称、传递三大关系 \[ \begin{align} a \equiv a \pmod m\\ a \equiv b \pmod m \iff b \equiv a \pmod m\\ a \equiv b \pmod m \and b \equiv c \pmod m \Rightarrow a \equiv c \pmod m \end{align} \] 加、减、乘与一般运算一致. 若\(a \equiv b \pmod m\)且\(c \equiv d \pmod m\), 则有 \[ a+c \equiv b+d \pmod m\\ a-c \equiv b-d \pmod m\\ a*c \equiv b*d \pmod m \] 除法比较特殊. 若\(ac \equiv bc \pmod m\), 则 \[ a \equiv b \pmod{\frac{m}{\gcd(c,m)}} \]
其它性质:
\[ \begin{align} ax \equiv 0 \pmod p \iff x \equiv 0 \pmod{\frac{p}{\gcd(p,a)}} \end{align}\\ x \% b == 0 \and \frac{x}{b} \equiv 0 \pmod p \iff x \equiv 0 \pmod{bm}\\ a \equiv b \pmod m \iff a^n \equiv b^n \pmod m\\ a \equiv b \pmod m \and n|m \Rightarrow a \equiv b \pmod n \]
## 快速幂
计算\(x^d \mod p\)
注意\(a^b % p = ((a % p)^b) % p\) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
ll Pow(ll x,ll d){
	ll tans = 1;
	if(d == 0)return 1%p;
    x %= p;
	ll a = Pow(x,d/2);
	tans = a*a%p;
	if(d%2)tans = tans*x%p; 
	return tans%p;
}

LL gcd(LL a,LL b)
{
    while(b)
    {
        LL tmp=b;
        b=a%b;
        a=tmp;
    }
    return a;
}
## 一次不定方程
\(\sum_{i = 1}^n a_i \cdot x_i = c\)有整数解当且仅当\(\gcd(a_1,a_2,...,a_n) | c\)
## GCD & exGCD
x,y,d为全局变量.求x和y使得 \(ax + by = d\) 且 \(|x| + |y|\)最小.
其中 \(d = gcd(a,b)\) 
1
2
3
4
5
6
7
8
9
10
11
ll gcd(ll a,ll b){
	return b?gcd(b,a%b):a;
}
void gcd(ll a,ll b,ll & d,ll & x,ll & y){
	//notice d = gcd(a,b)
	if(!b)d = a,x = 1,y = 0;
	else{
		gcd(b,a%b,d,y,x);
		y -= x*(a/b);
	}
}
## 逆元
返回模p下a的逆.不存在则返回-1. 
1
2
3
4
5
ll inv(ll a,ll p){
	ll d,x,y;
	gcd(a,p,d,x,y);
	return d == 1?(x+p)%p:-1;
}
### 线性筛逆元:
返回模p下1...n的逆.保存在invs里 
1
2
3
4
5
6
void get_inv(){
	inv[1] = 1;
	for(int i = 2;i < maxn;i++){
		inv[i] = (ll)(p - p/i)*inv[p%i]%p;
	}
}
## a,b不互质时
\[ \frac{a}{b} \mod m = \frac{a \mod (mb)}{b} \]
证: \[ \begin{align} &\frac{a}{b} = km + x (x < m)\\ \Rightarrow \space &a =kbm + bx\\ \Rightarrow \space &a \mod (bm) = bx\\ \Rightarrow \space &\frac{a \mod (bm)}{b} = x\\ \Rightarrow \space &\frac{a \mod (bm)}{b} = \frac{a}{b} \mod m \end{align} \]
## 中国剩余定理(CRT)
求同余方程组 \[ \begin{cases} x &\equiv a_1 \pmod {m_1} \\ x &\equiv a_2 \pmod {m_2} \\ &\vdots \\ x &\equiv a_k \pmod {m_k} \\ \end{cases} \] 在\([0,\prod_{i = 1}^nm_i])\)中的解. 要求任意\(m_i\)两两互质
输入: int n 方程组个数, ll a[] ll m[],方程组的参数和模数
复杂度:\(O(nlogM)\), 其中\(M\)与\(m_i\)同阶
输出: 方程组在\([0,\prod_{i = 1}^nm_i])\)中的解.
前置: 逆元
调用 CRT(n,a,m) 
1
2
3
4
5
6
7
8
9
10
11
ll CRT(int n,ll * a,ll * m){
	ll M = 1;
	for(int i = 1;i <= n;i++)M *= m[i];
	ll ans = 0;
	for(int i = 1;i <= n;i++){
		ll mt = M/m[i];
		ll ret = inv(mt,m[i]);
		ans = (ans + a[i]*mt%M*ret%M)%M;
	}
	return (ans + M)%M;
}
## Eratosthenes筛 
1
2
3
4
5
6
7
8
void prime(){
	for(int i = 2;i <= n;i++){
		if(!n_prime[i]){
			for(int j = i*2;j <= n;j += i)n_prime[j] = 1;
		}
	}
	n_prime[1] = 1;
}
## 线性筛
筛素因子时记得直接去掉最后一个素因子, 否则复杂度不对
### 线性筛素数 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
bool not_p[maxn];
vector<int>prime;
void euler_sieve(int n){
	not_p[1] = 1;
	for(int i = 2;i < n;i++){
		if(!not_p[i]){
			prime.push_back(i);
		}
		for(int p:prime){
			if((ll)p*i >= n)break;
			not_p[p*i] = 1;
			if(i%p == 0)break;
		}
	}
}
### 线性筛\(\varphi\)
\(\varphi(n)\): 小于等于\(n\)的正整数中与\(n\)互质的数个数
\(\varphi\)的一些性质(以下\(p\)代表素数) \[ \varphi(p) = p-1\\ 如果i \% p == 0, \varphi(i*p) = \varphi(i)*p\\ 否则\varphi(i*p) = \varphi(i)*(p-1) \]
phi[n]为\(\varphi(n)\),prime里存的素数 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
int phi[maxn];
bool not_p[maxn];
vector<int>prime;
void eular_sieve(int n){
	not_p[1] = phi[1] = 1;
	for(int i = 2;i < n;i++){
		if(!not_p[i]){
			phi[i] = i-1;
			prime.push_back(i);
		}
		for(int p:prime){
			if((ll)p*i >= n)break;
			not_p[p*i] = 1;
			if(i%p)phi[i*p] = phi[i]*(p-1);
			else{
				phi[i*p] = phi[i]*p;
				break;
			}
		}
	}
}
### 线性筛\(\mu\)
\[ \mu(n)= \begin{cases} 1&n=1\\ 0&n\text{ 含有平方因子}\\ (-1)^k&k\text{ 为 }n\text{ 的本质不同质因子个数}\\ \end{cases} \] 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
int phi[maxn],mu[maxn];
bool not_p[maxn];
vector<int>prime;
void sieve_mu(){
	not_p[1] = mu[1] = 1;
	for(int i = 2;i < maxn;i++){
		if(!not_p[i]){
			mu[i] = -1;
			prime.push_back(i);
		}
		for(int p:prime){
			if(ll(p)*i >= maxn)break;
			not_p[p*i] = 1;
			if(i%p)mu[i*p] = -mu[i];
			else{
				mu[i*p] = 0;
				break;
			}
		}
	}
}
## 杜教筛
详细请见杜教筛笔记.
杜教筛快速筛\(\varphi\)和\(\mu\)的前缀和.\(maxn = n^{\frac{2}{3}}\)时复杂度为\(O(n^{\frac{2}{3}})\) 
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
typedef long long ll;
const int maxn = 1664510;
bool not_p[maxn];
vector<int>prime;
ll mu[maxn],phi[maxn];
void sieve(){
	phi[1] = mu[1] = 1;
	for(int i = 2;i < maxn;i++){
		if(!not_p[i]){
			mu[i] = -1;
			phi[i] = i-1;
			prime.push_back(i);
		}
		for(int p:prime){
			if(1LL*p*i >= maxn)break;
			not_p[i*p] = 1;
			if(i%p){
				phi[i*p] = phi[i]*(p-1);
				mu[i*p] = -mu[i];
			}
			else{
				phi[i*p] = phi[i]*p;
				mu[i*p] = 0;
				break;
			}
 		}
	}
	for(int i = 1;i < maxn;i++)mu[i] += mu[i-1],phi[i] += phi[i-1];
}
map<int,ll>Mu,Phi;
ll du_mu(int x){
	if(x < maxn)return mu[x];
	else if(Mu[x])return Mu[x];
	ll tans = 1;
	for(int l = 2,r = 0;l <= x;l = r+1){
		r = x/(x/l);
		tans -= 1LL*(r - l + 1)*du_mu(x/l);
	}
	return Mu[x] = tans;
}
ll du_phi(int x){
	if(x < maxn)return phi[x];
	else if(Phi[x])return Phi[x];
	ll tans = 1LL*x*(x+1)/2;
	for(int l = 2,r = 0;l <= x;l = r+1){
		r = x/(x/l);
		tans -= 1LL*(r-l+1)*du_phi(x/l);
	}
	return Phi[x] = tans;
}
## 欧拉函数\(φ(x)\) 
1
2
3
4
5
6
7
8
9
10
ll cal_phi(ll x){
	int m = (ll)sqrt(x+0.5);
	ll tans = x;
	for(int i = 2;i <= m;i++)if(!(x%i)){
		tans = tans/i*(i-1);
		while(!(x%i))x /= i;
	}
	if(x > 1)tans = tans / x * (x-1);
	return tans;
}
### 线性筛\(φ(x)\) 
1
2
3
4
5
6
7
8
void sieve_phi(int n){
	phi[1] = 1;
	for(int i = 2;i <= n;i++)if(!phi[i])
		for(int j = i;j <= n;j += i){
			if(!phi[j])phi[j] = j;
			phi[j] = phi[j]/i * (i-1);
		}
}
## 有理数取模
求\(a/b\) % \(p\) （p为质数)
由费马小定理\(b^{p-1} ≡ 1 \pmod p\) 可得 \(b^{p-2} ≡ b^{-1} \pmod p\)
故只需求\(a*b^{p-2} \bmod p\)
a,b过大则先膜一下
## 欧拉定理
若\(a,p\)互质,则有
\(a^{φ(p)} ≡ 1 \pmod p\)
### 扩展欧拉定理(无需a,p互质)
\(b ≥ φ(p)\)时:
\(a^b ≡ a^{b \space \bmod \space φ(p) +φ(p)} \pmod p\)
例题: 给出\(w_1,w_2,..w_n\)和\(q\)组询问\(l \space r\). 求\(w_l^{w_{l+1}^{...^{w_r}}} \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
//CF906D 
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll maxn = 1e5+10;
ll raw[maxn],n,mod;
long long r1aw[maxn];
map<ll,ll>vis;
ll phi(ll x){
	ll m = (ll)sqrt(x+0.5);
	ll tans = x;
	for(ll i = 2;i <= m;i++)if(!(x%i)){
		tans = tans/i*(i-1);
		while(!(x%i))x /= i;
	}
	if(x > 1)tans = tans / x * (x-1);
	return tans;
}
ll Pow(ll x,ll d,ll p){
	ll tans = 1;
	if(d == 0)return 1%p;
    x %= p;
	ll a = Pow(x,d/2,p);
	tans = a*a%p;
	if(d%2)tans = tans*x%p; 
	return tans%p;
}
ll w[maxn],p0[maxn];
int cal(int l,int r,int k) {
	if (l==r) return min(w[l],w[l]%p0[k]+p0[k]);
	if (p0[k]==1) {
		return 1;
	} else {
		int c=cal(l+1,r,k+1);
		if (w[l]==1) return 1;
		ll v=Pow(w[l],c,p0[k]);
		ll d=1ll<<60,e=1;
		for(int i = 0;i < c;i++) { d/=w[l]; e*=w[l]; if (d==0) break; }
		if (d==0) return v+p0[k]; else return min(e,e%p0[k]+p0[k]);
	}
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	int n;
	cin >> n >> mod;
	int idx = 0;
	p0[idx] = mod;
	while(p0[idx] != 1){p0[idx+1] = phi(p0[idx]);idx++;}
	for(ll i = 1;i <= n;i++){
		cin >> w[i];
	}
	int q;
	cin >> q;
	while(q--){
		int l,r;
		cin >> l >> r;
		cout << cal(l,r,0)%mod << endl;
	}
	return 0;
}
## 原根
proot(p) 返回素数p的原根 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int>pri;
bool g_test(ll g,ll p){
	for(ll x:pri){
		if(Pow(g,(p-1)/x,p) == 1)return 0;
	}
	return 1;
}
ll proot(ll p){
	ll tmp = p-1;
	for(ll i = 2;i <= tmp/i;i++){
		if(tmp%i == 0){
			pri.push_back(i);
			while(tmp%i == 0)tmp /= i;
		}
	}
	if(tmp != 1)pri.push_back(tmp);
	ll g = 1;
	while(1){
		if(g_test(g,p))return g;
		g++;
	}
}
## 离散对数(BSGS)
### 原版
只支持模数为\(p\)的情况.
bsgs(ll a,ll b,ll p)返回最小的\(x\)使得\(a^x \equiv b\pmod{p}\), 不存在则返回\(-1\). 要求\(p\)为素数. 时间复杂度\(O(\sqrt p \cdot \log(\sqrt{p}))\) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int bsgs(int a,int b,int p){
	map<int,int>f;
	int m = ceil(sqrt(p));
	b %= p;
	for(int i = 1;i <= m;i++){
		b = 1ll*b*a%p;
		f[b] = i;
	}
	int v = Pow(a,m,p);
	b = 1;
	for(int i = 1;i <= m;i++){
		b = 1ll*b*v%p;
		if(f.count(b))return (1ll*i*m - f[b] + p)%p;
	}
	return -1;
}
求最小的正整数\(n\)使得\(a^n \equiv 1 \pmod m\) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
int ord(int a,int m){
	a %= m;
	if(__gcd(a,m) != 1)return -1;
	int sq = sqrt(m)+1;
	map<int,int>mp;
	ll s = a;
	for(int i = 1;i <= sq;i++){
		if(!mp.count(s))mp[s] = i;
		s = s*a%m;
	}
	ll g = inv(Pow(a,sq,m),m),x = 1;
	for(int k = 0;k <= m/sq;k++){
		if(mp.count(x))return sq*k + mp[x];
		x = x*g%m;
	}
}
### 扩展版
\(p\)可以是任意模数. 前置:bsgs exgcd inv 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int exbsgs(int a,int b,int p){
	if(b == 1 or p == 1)return 0;
	int g = __gcd(a,p),k = 0,na = 1;
	while(g > 1){
		if(b%g != 0)return -1;
		k++;
		b /= g;p /= g;na = na*(a/g)%p;
		if(na == b)return k;
		g = __gcd(a,p);
	}
	int f = bsgs(a,b*inv(na,p)%p,p);
	if(f == -1)return -1;
	return f+k;
}
## 拉格朗日插值
px,py为点坐标,p为模数,k是L(k)那个k inv为逆元
返回L(k) 
1
2
3
4
5
6
7
8
9
10
11
12
13
ll lagrange(ll n,ll k){
	ll tans = 0;
	for(int i = 1;i <= n;i++){
		ll f = 1,g = 1;
		for(int j = 1;j <= n;j++)if(i != j){
			f = f*(k - px[j]+p)%p;
			g = g*(px[i] - px[j]+p)%p;
		}
		ll lota = py[i]*f%p*inv(g,p)%p;
		tans = (tans+lota)%p;
	}
	return tans; 
}
### 线性求\(\sum_{i = 0}^{n}i^k\)
(返回L(k)的值,多项式次数为n.按需更改yi的值,复杂度\(O(k)\))
(注意求\(\sum_{i = 0}^{n}i^k\)需要调用line_lagrange(k+2,n) ) 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
ll line_lagrange(ll n,ll k){
	ll tans = 0;
	get_inv();
	pre[0] = suf[n+1] = ifac[0] =  1;
	for(int i = 1;i <= n;i++){
		pre[i] = pre[i-1]*(k-i)%mod;
		ifac[i] = ifac[i-1]*invs[i]%mod;
	}
	for(int i = n;i >= 0;i--)suf[i] = suf[i+1]*(k-i)%mod;
	ll yi = 0; 
	for(int i = 1;i <= n;i++){
		yi = (yi+pow_mod(i,n-2))%mod;
		ll a = pre[i-1]*suf[i+1]%mod*yi%mod*ifac[i-1]%mod*ifac[n-i]%mod;
		if((n-i)%2)a = mod-a; 
		tans = (tans+a+mod)%mod;
	}
	return tans;
}
## 卢卡斯定理
Lucas(n,m) = C(n%p,m%p)*Lucas(n/p,m/p)%p
p需为素数,fac[]为预处理的阶乘,时间复杂度大约是\(O(p + mlogm)\)(预处理需要\(O(p)\)) 
1
2
3
4
5
6
7
8
ll C(ll n,ll m){
	if(m > n)return 0;
	return (fac[n]*Pow(fac[m],p-2)%p*Pow(fac[n-m],p-2)%p)%p;
}
ll Lucas(ll n,ll m){
	if(!m)return 1;
	return C(n%p,m%p)*Lucas(n/p,m/p)%p;
}

线性基

maxbit取值域位数

insert返回0则说明现有基可以表示x

qmax返回异或最大值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
struct LinerBasis{
	ll base[maxbit+4];
	bool vis[maxbit+4];
	bool insert(ll x){
		for(int i = maxbit;~i;i--){
			if(x&(1LL<<i)){
				if(!base[i]){
					base[i] = x;
					break;
				}
				else x ^= base[i];
			}
		}
		return x;
	}
	ll qmax(){
		ll ans = 0;
		for(int i = maxbit;~i;i--)ans = max(ans,ans^base[i]);
		return ans;
	}
}yay;

杜教BM(快速求模意义下线性递推式)

设\(n\)为输入项数, 复杂度\(O(n^2\log n)\). 对于\(k\)阶递推式(例如斐波那契数列为2阶递推式)至少需要前\(2k\)项

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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include<bits/stdc++.h>
using namespace std;
#define rep(i,a,n) for (int i=a;i<n;i++)
#define pb push_back
typedef long long ll;
#define SZ(x) ((ll)(x).size())
typedef vector<ll> VI;
typedef pair<ll, ll> PII;
const ll mod = 1000000007;
ll powmod(ll a, ll b) { ll res = 1; a %= mod; assert(b >= 0); for (; b; b >>= 1) { if (b & 1)res = res * a % mod; a = a * a % mod; }return res; }
// head
 
ll _, n;
namespace linear_seq {
    const ll N = 10010;
    ll res[N], base[N], _c[N], _md[N];
 
    vector<ll> Md;
    void mul(ll* a, ll* b, ll k) {
        rep(i, 0, k + k) _c[i] = 0;
        rep(i, 0, k) if (a[i]) rep(j, 0, k) _c[i + j] = (_c[i + j] + a[i] * b[j]) % mod;
        for (ll i = k + k - 1; i >= k; i--) if (_c[i])
            rep(j, 0, SZ(Md)) _c[i - k + Md[j]] = (_c[i - k + Md[j]] - _c[i] * _md[Md[j]]) % mod;
        rep(i, 0, k) a[i] = _c[i];
    }
    ll solve(ll n, VI a, VI b) { // a 系数 b 初值 b[n+1]=a[0]*b[n]+...
//        printf("%d\n",SZ(b));
        ll ans = 0, pnt = 0;
        ll k = SZ(a);
        assert(SZ(a) == SZ(b));
        rep(i, 0, k) _md[k - 1 - i] = -a[i]; _md[k] = 1;
        Md.clear();
        rep(i, 0, k) if (_md[i] != 0) Md.push_back(i);
        rep(i, 0, k) res[i] = base[i] = 0;
        res[0] = 1;
        while ((1ll << pnt) <= n) pnt++;
        for (ll p = pnt; p >= 0; p--) {
            mul(res, res, k);
            if ((n >> p) & 1) {
                for (ll i = k - 1; i >= 0; i--) res[i + 1] = res[i]; res[0] = 0;
                rep(j, 0, SZ(Md)) res[Md[j]] = (res[Md[j]] - res[k] * _md[Md[j]]) % mod;
            }
        }
        rep(i, 0, k) ans = (ans + res[i] * b[i]) % mod;
        if (ans < 0) ans += mod;
        return ans;
    }
    VI BM(VI s) {
        VI C(1, 1), B(1, 1);
        ll L = 0, m = 1, b = 1;
        rep(n, 0, SZ(s)) {
            ll d = 0;
            rep(i, 0, L + 1) d = (d + (ll)C[i] * s[n - i]) % mod;
            if (d == 0) ++m;
            else if (2 * L <= n) {
                VI T = C;
                ll c = mod - d * powmod(b, mod - 2) % mod;
                while (SZ(C) < SZ(B) + m) C.pb(0);
                rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod;
                L = n + 1 - L; B = T; b = d; m = 1;
            }
            else {
                ll c = mod - d * powmod(b, mod - 2) % mod;
                while (SZ(C) < SZ(B) + m) C.pb(0);
                rep(i, 0, SZ(B)) C[i + m] = (C[i + m] + c * B[i]) % mod;
                ++m;
            }
        }
        return C;
    }
    ll gao(VI a, ll n) {
        VI c = BM(a);
        c.erase(c.begin());
        rep(i, 0, SZ(c)) c[i] = (mod - c[i]) % mod;
        return solve(n, c, VI(a.begin(), a.begin() + SZ(c)));
    }
};
int main() {
    while (~scanf("%lld", &n)) {
        vector<ll>v;
        v.push_back(1);/*前几项*/v.push_back(5);v.push_back(15);
        v.push_back(35);v.push_back(70);v.push_back(126);
        v.push_back(210);v.push_back(330);v.push_back(495);v.push_back(715);
        v.push_back(1001);
        //输入n ,输出第n项的值  一般大于10项即可出答案，越多越好
        printf("%lld\n", linear_seq::gao(v, n - 1));
    }
}

特殊计数序列

第二类Stirling数

第二类Stirling数\(S(n,m)\)计数的是把\(n\)个元素的集合划分到\(m\)个不可区分的盒子中且没有空盒子的划分个数.

\(S(n,n) = 1\), \(S(n,0) = 0\space(n > 0)\)

递推式: \(S(n,m) = m \cdot S(n-1,m) + S(n-1,m-1)\)

求第n行

\[ S(n,m) = \sum_{i = 0}^m \frac{i^n}{i!}\cdot\frac{(-1)^{m-i}}{(m-i)!} \] 用NTT求解即可. finv[i]为\(i!\)的逆元, 需要线性预处理出来. mul()见NTT板子. 时间复杂度\(O(n \log n)\)

1
2
3
4
5
6
7
8
9
10
11
Poly Stirling_row(int n){
	Poly A(n+1),B(n+1);
	for(int i = 0;i <= n;i++){
		A[i] = 1ll*Pow(i,n)*finv[i]%p;
		B[i] = ( ( 1ll * ((n-i)%2?-1:1) * finv[n-i] )%p+p)%p;
	}
	reverse(B.begin(),B.end());
	auto tans = mul(A,B);
	tans.resize(n+1);
	return tans;
}

求第k列

\[ \sum_{n = k}^{\infin}S(n,k) \cdot \frac{x^n}{n!} = \frac{(e^x - 1)^k}{k!} \]

因而有 \[ S(n,k) = \frac{n!}{k!}\cdot[n](e^x - 1)^k \] 多项式快速幂即可. 时间复杂度\(O(n \log n)\)

1
2
3
4
5
6
7
8
Poly Stirling_Col(int n,int k){
	Poly A(n+1);
	for(int i = 1;i <= n;i++)A[i] = finv[i];
	A = powP(A,k,k);
	A.resize(n+1);
	for(int i = 0;i <= n;i++)A[i] = 1ll*A[i]*finv[k]%p*fac[i]%p;
		return A;
}

Bell数

Bell数(\(B(n)\))计数的是将一个有\(n\)个元素的集合划分成任意个非空子集的方案数. 显然有\(B(n) = \sum_{i = 1}^n S(n,i)\)

前10个Bell数为1 2 5 15 52 203 877 4140 21147 115975 \[ B(n) = [n]e^{e^x - 1} \] 多项式指数函数即可. 时间复杂度\(O(n\log n)\)

1
2
3
4
5
6
7
Poly Bell(int n){
	Poly A(n+1);
	for(int i = 1;i <= n;i++)A[i] = finv[i];
	A = expP(A);
	for(int i = 0;i <= n;i++)A[i] = 1ll*A[i]*fac[i]%p;
	return A;
}

博弈论

威佐夫游戏

两人轮流取两堆筹码，其中取法有两种：取走一堆中任意个筹码，或从两堆中取走相同数目的筹码。取完所有筹码的一方获胜。

(1为先手胜)

1
2
3
4
if(n1>n2)  swap(n1,n2);
temp=floor((n2-n1)*(1+sqrt(5.0))/2.0);
if(temp==n1) cout<<"0"<<endl;
else cout<<"1"<<endl;

矩阵

定义&快速幂

注意：vector版的矩阵会比数组慢起码2倍，时间吃紧的时候慎用

数组版:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#define FOR(i,a,b) for(int i = a;i <= b;i++)
struct Matrix{
	ll A[maxn][maxn];
	void initer(){
        for(int i = 1;i <= n;i++)
            for(int j = 1;j <= n;j++)A[i][j] = 0;
		for(int i = 1;i <= n;i++)A[i][i] = 1;
	}
	Matrix operator * (const Matrix & a){
		Matrix res;
		memset(res.A,0,sizeof(res.A));
		FOR(k,1,n)
			FOR(i,1,n)
				FOR(j,1,n){
					res.A[i][j] = (res.A[i][j] + (A[i][k]*a.A[k][j]));
					if(res.A[i][j] > MOD2)res.A[i][j] -= MOD2;
				}
		FOR(i,1,n)
			FOR(j,1,n)res.A[i][j] %= MOD;
		return res;
	}
}

vector版:

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
typedef vector<vector<int> > Matrix;
void Minit(Matrix & a){
	int n = a.size();
	for(int i = 0;i < n;i++){
		for(int j = 0;j < n;j++)a[i][j] = (i == j);
	}
}
Matrix operator * (const Matrix & a,const Matrix & b){//需保证a的列数等于b的行数且两个矩阵都不为空
	int n = a.size(),m = b.size(),r = b[0].size();
	Matrix res(n,vector<int>(m));
	for(int k = 0;k < m;k++){
		for(int i = 0;i < n;i++){
			for(int j = 0;j < r;j++){
				res[i][j] = (res[i][j] + (1ll*a[i][k]*b[k][j])%MOD)%MOD;
			}
		}
	}
	return res;
}
Matrix matrixPow(Matrix x,ll d){
	int n = x.size();
	Matrix ans(n,vector<int>(n));
	Minit(ans);
	while(d){
		if(d&1)ans = ans*x;
		x = x*x;
		d /= 2;
	}
	return ans;
}

置换矩阵

每一行每一列有且仅有一个\(1\), 其余元素为\(0\)的矩阵为置换矩阵.

置换矩阵的逆是它的转置.

循环矩阵

除第一行外, 其余每行都是上一行循环右移一格的矩阵是循环矩阵. 循环矩阵的乘积仍是循环矩阵.

因此只需保存和计算循环矩阵的第一行, 两个循环矩阵相乘复杂度为\(O(n^2)\)

image-20210924200136207 
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
typedef vector<int> cirMatrix;//循环矩阵只需保存第一行
Matrix buildCirMatrix(const cirMatrix & a){//构建循环矩阵,此处为右移
	int n = a.size();
	Matrix res(n,vector<int>(n));
	res[0] = a;
	for(int i = 1;i < n;i++){
		res[i][0] = res[i-1][n-1];
		for(int j = 1;j < n;j++)res[i][j] = res[i-1][j-1];
	}
	return res;
}
cirMatrix operator * (const cirMatrix & a,const cirMatrix & b){
	Matrix A = Matrix(1,a);
	Matrix B = buildCirMatrix(b);
	return (A*B)[0];
}
cirMatrix matrixPow(cirMatrix x,ll d){
	int n = x.size();
	cirMatrix ans(n);
	ans[0] = 1;
	while(d){
		if(d&1)ans = ans*x;
		x = x*x;
		d /= 2;
	}
	return ans;
}

高斯消元

输入一个\(n\)行\(n+1\)列的矩阵,其中\(A[i][n]\)代表等式右边第\(i\)个值. 复杂度\(O(n^3)\)

gauss(矩阵,n)返回1有解, 返回0无解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
bool gauss(vector<vector<ld> > & A,int n){
	for(int i = 0;i < n;i++){
		int r = i;
		for(int j = i+1;j < n;j++)if(fabs(A[j][i]) > fabs(A[r][i]))r = j;
		if(r != i)for(int j = 0;j < n+1;j++)swap(A[i][j],A[r][j]);
		if(fabs(A[i][i]) <= eps)return 0;
		for(int j = n;j >= i;j--){
			for(int k = i+1;k < n;k++){//提高精度
				A[k][j] -= A[k][i]/A[i][i]*A[i][j];
			}
		}
	}
	for(int i = n-1;i >= 0;i--){
		for(int j = i+1;j < n;j++)A[i][n] -= A[j][n]*A[i][j];
		A[i][n] /= A[i][i];
	}
	return 1;
}

异或方程组

时间复杂度\(O(\frac{n^3}{w})\), \(w = 32或64\), 取决于评测机位数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
const int maxb = 1005;
vector<bitset<maxb> > a(m,bitset<maxb>());
int xorRank(vector<bitset<maxb> > & a,int n,int m){
	//n个未知数, m个方程
	int mr = -1,i = 0,j = 0;
	while(i < m and j < n){
		int r = i;
		for(int k = i;k < m;k++)if(a[k][j]){r = k;break;}
		if(a[r].test(j)){
			mr = max(mr,r);
			if(r != i)swap(a[r],a[i]);
			for(int u = i+1;u < m;u++)if(a[u].test(j))a[u] ^= a[i];
			i++;
		}
		//else return 0; //某一行全为0, 说明有多解
		j++;
	}
	for(int k = m-1;k >= 0;k--){
		for(int j = k+1;j < m;j++)a[k][m-1] = (a[k][m-1] ^ (a[j][m-1] & a[k][j]));
	}
    return i;//返回该矩阵的秩
	//return mr+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
const ld eps = 1e-9;
ld caldet(vector<vector<ld> > & a,int n){
	ld det = 1;
	for(int i = 0;i < n;i++){
		int k = i;
		for(int j = i+1;j < n;j++){
			if(fabs(a[j][i]) > fabs(a[k][i]))k = j;
		}
		if(fabs(a[k][i]) < eps){
			det = 0;
			break;
		}
		swap(a[i],a[k]);
		if(i != k)det = -det;
		det *= a[i][i];
		for(int j = i+1;j < n;j++)a[i][j] /= a[i][i];
		for(int j = 0;j < n;j++){
			if(j != i and fabs(a[j][i]) > eps){
				for(int k = i+1;k < n;k++){
					a[j][k] -= a[i][k] * a[j][i];
				}
			}
		}
	}
	return det;
}

整数版

输入一个vector矩阵, 输出行列式模\(p\)的值. 均摊复杂度\(O(n^3 + n^2\log n)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int caldet(vector<vector<int> > & a,int n,int p){
	int det = 1,w = 1;
	for(int i = 0;i < n;i++){
		for(int j = i+1;j < n;j++){
			while(a[i][i]){
				int div = a[j][i]/a[i][i];
				for(int k = i;k < n;k++){
					a[j][k] = (a[j][k] - 1ll*div*a[i][k]%p + p)%p;
				}
				swap(a[i],a[j]);w = -w;
			}
			swap(a[i],a[j]);w = -w;
		}
	}
	for(int i = 0;i < n;i++)det = 1ll*a[i][i]*det % p;
	det = (1ll*w*det + p)%p;
	return det;
}

Kirchhoff 矩阵树定理

矩阵树定理用于求解图的生成树个数. 来源: https://oi-wiki.org/graph/matrix-tree/

设\(t(G)\)为图\(G\)的全部生成树边权积的和(当边权为1时即为全部生成树的个数)

即 \[ t(G) = \sum_{T \in G} \prod_{e \in T}w_e \] 其中\(T\)代表图\(G\)的生成树, \(e\)为生成树上的边，\(w_e\)为边的权值.

以下所讨论的图允许重边但不允许自环(去掉自环即可).

无向图形式-行列式

设\(D\)为无向图\(G\)的度数矩阵, \(A\)为无向图的邻接矩阵. \(deg(x)\)为点\(x\)的度数.

则\(D[i][i] = deg(i) \cdot[i == j]\), \(A[i][j] = |i与j间的连边权值之和|\)

设\(Laplace\)矩阵\(L = D - A\). \(L_r\)为\(L\)的\(n-1\)阶主子式(即将矩阵\(L\)去掉第\(r\)行第\(r\))列.

则对任意的\(r\), \(L_r\)的行列式为该无向图的生成树个数. 即\(t(G) = det(L_r)\)

无向图形式-特征值

若\(\lambda_1,\lambda_2,...,\lambda_{n-1}\)为\(L\)的\(n-1\)个非负特征值, 则有 \[ t(G) = \frac{1}{n}\prod_{i = 1}^{n-1}\lambda_i \]

有向图形式

有向图的生成树分两种: 根向, 所有边均由儿子指向父亲; 叶向, 所有边均由父亲指向儿子.

设\(D^{out}\)为有向图\(G\)的出度矩阵, \(D^{in}\)为有向图\(G\)的入度矩阵, 两者与无向图形式相同.

\(A\)为有向图邻接矩阵, \(A[i][j] = |i指向j的边数|\)

设入度Laplace矩阵\(L^{in} = D^{in} - A\), 出度Laplace矩阵为\(L^{out} = D^{out} - A\)

则以\(k\)为根的叶向生成树个数为\(t^{leaf}(G,k) = det(L_k^{in})\), 根向生成树个数为\(t^{root}(G,k) =det(L_k^{out})\).

若要求所有生成树, 枚举\(k\)并求和即可.

带权形式

对有向图无向图均适用. 将带权边转化为无权的重边, 重边次数为边权即可.

BEST定理

通过图中所有边恰好一次且行遍所有顶点的回路称为欧拉回路.

具有欧拉回路的无向图或有向图称为欧拉图.

设 \(G\) 是有向欧拉图，那么 \(G\) 的不同欧拉回路总数 \(ec(G)\) 是 \[ ec(G) = t^{root}(G,k)\prod_{v\in V}(\deg (v) - 1)! \]

注意，对欧拉图 \(G\) 的任意两个节点 \(k, k'\)，都有 \(t^{root}(G,k)=t^{root}(G,k')\)，且欧拉图 \(G\) 的所有节点的入度和出度相等。

一些例题

加速递推式

以\(Fibonacci\)为例:

\(f(1) = f(2) = 1,f(n) = f(n-1) + f(n-2),(n \geq 10^9)\),求\(f(n)\)

设 \[ F_i = \begin{bmatrix} f_{i+1}\\ f_i \end{bmatrix}, F_0 = \begin{bmatrix} f_1\\ f_0 \end{bmatrix} \] 且\(F_i = A^i \cdot F_0\)

则有 \[ A = \begin{bmatrix} 1 & 1\\ 1 & 0 \end{bmatrix} \] 一般地, 对于递推式\(f(i) = \sum_{i = 1}^{k}a_i \cdot f(i-k)\)

有: \[ F_i = \begin{bmatrix} f(i+k-1)\\ f(i+k-2)\\ ...\\ f(i) \end{bmatrix}, A = \begin{bmatrix} a_1 & a_2 & a_3 & ... & a_k\\ 1 & 0 & 0& \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ 0 & 0 & 1 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix} \] 且\(F_n = A^n \cdot F_0\)

求包含严格k条边的方案数/最短路

1.无权图求从\(i\)到\(j\)经过\(k\)条边的方案数:求邻接矩阵的\(k\)次幂即可

2.有向图求包含\(k\)条边的最短路:求邻接矩阵的\(k\)次幂,把普通矩阵乘法改为Floyd

0-1分数规划

给出数列\(\{a_i\},\{b_i\}\),求一组\(w_i \in \{0,1\}\), 使得 \[ \frac{\sum_{i = 1}^na_i \cdot w_i}{\sum_{i = 1}^nb_i \cdot w_i} \] 最大或最小

或给出\(n\)个物品,每个物品有两个权值\(a_i\)和\(b_i\), 选择任意个物品使得\(\frac{\sum a_i}{\sum b_i}\)最大

二分法

\[ \frac{\sum_{i = 1}^na_i \cdot w_i}{\sum_{i = 1}^nb_i \cdot w_i} \ge M \Leftrightarrow \sum_{i = 1}^n w_i \cdot (a_i - M\cdot b_i) \ge 0 \]

因此只需判断$_{i = 1}^n w_i (a_i - Mb_i) \(的最大值是否非负即可对\)M$二分.

Dinkelbach 算法(迭代法)

将上一轮的答案作为输入进行迭代.注意此时答案应为\(\frac{\sum_{i = 1}^na_i}{\sum_{i = 1}^nb_i}\), 而不是变形后的式子.

例题

[POJ2728]Desert King

给出\(n \leq 10^3\)个点的横纵坐标和高度, 两点间距离为欧几里得距离, 价值为高度之差绝对值. 所有点构成一个完全图. 求一棵生成树使得

距离/价值比最小.

二分法

此题是有一定限制的0-1分数规划(所选的"物品"构成一棵生成树). 直接求边权为\(w_i \cdot (a_i - M\cdot b_i)\)的最小生成树, 判断是否有解即可二分.

由于是完全图, 采用不加优化的Prim求解.

复杂度大概是\(O(n^2 \log((R - L)*2^{-eps}))\)

由于二分法较慢, 本题需要猜一个答案上界(100)才能通过=-=

用时2610ms

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
//POJ2728
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
typedef double ld;
#define endl '\n'
const int maxn = 1e3+10;
const ld INF = 1e6+10;
const ld eps = 1e-5;
ld G[maxn][maxn],dis[maxn];
int x[maxn],y[maxn],h[maxn],n;
ld w[maxn][maxn];
bool vis[maxn];
ld Dis(int idx1,int idx2){
	return sqrt((x[idx1] - x[idx2])*(x[idx1] - x[idx2])+(y[idx1] - y[idx2])*(y[idx1] - y[idx2]));
}
bool prim(ld M){
	for(int i = 0;i <= n;i++)vis[i] = 0,dis[i] = INF;
	for(int i = 1;i < n;i++){
		int x = 0;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(!x or dis[x] > dis[j])x = j;
		}
		vis[x] = 1;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			dis[j] = min(dis[j],w[j][x] - M*G[j][x]);
		}
	}
	ld tans = 0;
	for(int i = 2;i <= n;i++)tans += dis[i];
	return tans >= 0.0;
}
signed main(){
	while(scanf("%d",&n) == 1){
		if(!n)break;
		for(int i = 1;i <= n;i++){
			scanf("%d%d%d",&x[i],&y[i],&h[i]);
		}
		for(int i = 1;i <= n;i++){
			for(int j = 1;j <= n;j++){
				G[i][j] = Dis(i,j);
				w[i][j] = abs(h[i] - h[j]);
			}
		}
		ld M = 0,L = 0,R = 100;
		while(R - L >= eps){
			M = (L+R)/2;
			if(prim(M))L = M;
			else R = M;
		}
		printf("%.3f\n",L);
	}
	return 0;
}

迭代法

求生成树时仍按边权\(w_i \cdot (a_i - M\cdot b_i)\)来计算. 但统计答案则计算\(\frac{\sum_{i = 1}^na_i}{\sum_{i = 1}^nb_i}\)

跑的飞快(250ms)

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
//POJ2728
#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long ll;
typedef double ld;
#define endl '\n'
const int maxn = 1e3+10;
const ld INF = 1e9+10;
const ld eps = 1e-5;
ld G[maxn][maxn],dis[maxn];
int x[maxn],y[maxn],h[maxn],to[maxn],n;
ld w[maxn][maxn];
bool vis[maxn];
ld Dis(int idx1,int idx2){
	return sqrt((x[idx1] - x[idx2])*(x[idx1] - x[idx2])+(y[idx1] - y[idx2])*(y[idx1] - y[idx2]));
}
ld prim(ld M){
	ld up = 0,down = 0;
	for(int i = 0;i <= n;i++)vis[i] = 0,dis[i] = INF,to[i] = 1;
	for(int i = 1;i < n;i++){
		int x = 0;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(!x or dis[x] > dis[j])x = j;
		}
		vis[x] = 1;
		for(int j = 1;j <= n;j++)if(!vis[j]){
			if(dis[j] > w[x][j] - M*G[x][j])dis[j] = w[x][j] - M*G[x][j],to[j] = x;
		}
	}
	for(int i = 2;i <= n;i++){
		up += w[i][to[i]];
		down += G[i][to[i]];
	}
	return up/down;
}
signed main(){
	while(scanf("%d",&n) == 1){
		if(!n)break;
		for(int i = 1;i <= n;i++){
			scanf("%d%d%d",&x[i],&y[i],&h[i]);
		}
		for(int i = 1;i <= n;i++){
			for(int j = 1;j <= n;j++){
				G[i][j] = Dis(i,j);
				w[i][j] = abs(h[i] - h[j]);
			}
		}
		ld ans = 0,lans = 1;
		while(fabs(lans - ans) > eps){
			lans = ans;
			ans = prim(ans);
		}
		printf("%.3f\n",ans);
	}
	return 0;
}

多项式

FFT

复数类操作:

1
2
3
complex<double>m;
m.real();m.imag();//获取实部虚部的值
m.real(114);m.imag(514);//为实部虚部赋值

包括手写复数类, DFT与IDFT, 多项式乘法.

  • 复数类: \(x\)为实部, \(y\)为虚部

  • FFT(Comp * A,int siz,int type)

    A: 复数数组, 全局变量, 为待求解的多项式

    siz: 多项式长度. 需保证为2的次方

    type: 1为DFT, -1为IDFT

  • Poly mul(Poly A,Poly B)

    多项式乘法. 输入两个多项式\(A,B\), 长度分别为\(n,m\)(也就是\(n-1,m-1\)阶). 返回一个长度\(n - m + 1\)的多项式, 为\(A,B\)的乘积.

    Poly为vector<int>.

!!!切记在点值表示法下相乘时, 两个IDFT的长度必须相同!!!

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
typedef double ld;
typedef vector<int> Poly;
const ld PI = acos(-1.0);
int rev[maxn];
struct Comp{
    double x, y;
    Comp (double x = 0, double y = 0):x(x),y(y){}
    Comp operator * (Comp B){
   		return Comp(x*B.x - y*B.y,x*B.y + y*B.x);
	}
	Comp operator + (Comp B){
	    return Comp(x + B.x,y + B.y);
	}
	Comp operator - (Comp B){
	    return Comp(x - B.x,y - B.y);
	}
};
Comp F1[maxn],F2[maxn];
void FFT(Comp * A,int siz,int type){
	int n = siz;
	int S = log2(n);
	for(int i = 0;i < n;i++)rev[i] = (rev[i>>1] >> 1) | ((i & 1) << (S - 1));
	for(int i = 0;i < n;i++)if(i < rev[i])swap(A[i],A[rev[i]]);
	for(int i = 1;i < n;i *= 2){
		Comp Wn(cos(PI/i),type*sin(PI/i));
		for(int j = 0;j < n;j += i*2){
			Comp W(1,0);
			for(int k = 0;k < i;k++){
				Comp facx = A[j+k],facy = W*A[j+k+i];
				A[j+k] = facx + facy;
				A[j+k+i] = facx - facy;
				W = W*Wn;
			}
		}
	}
    if(type == -1)for(int i = 0;i < n;i++)A[i].x = (int)((A[i].x/n + 0.5));
}
Poly mul(Poly A,Poly B){
	int n = A.size(),m = B.size();
	int siz = n + m - 1;
	Poly C(siz);
	if(siz < 64){//小于等于64项直接暴力算
		for(int i = 0;i < n;i++){
			for(int j = 0;j < m;j++)C[i+j] = (C[i+j] + 1LL*A[i]*B[j]);
		}
		return C;
	}
	int fsiz = 1;
	while(fsiz <= siz)fsiz *= 2;
	for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
	for(int i = 0;i < n;i++)F1[i] = A[i];
	for(int i = 0;i < m;i++)F2[i] = B[i];
	FFT(F1,fsiz,1);
	FFT(F2,fsiz,1);
	for(int i = 0;i < fsiz;i++)F1[i] = F1[i]*F2[i];
	FFT(F1,fsiz,-1);
	for(int i = 0;i < siz;i++){
		C[i] = ((ll)F1[i].x);
	}
	return C;
}

卷积

\[ C[i] = \sum_{j = 0}^iA[j]B[i-j] \]

卷积等同于多项式乘积

NTT

从OI-wiki抄来的, 用法和上面的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
const int maxn = 6e6+10;
const int p = 998244353;
int Pow(int x,int d){
	int tans = 1;
	if(d == 0)return 1%p;
	int a = Pow(x,d/2);
	tans = 1ll*a*a%p;
	if(d%2)tans = 1ll*tans*x%p; 
	return tans%p;
}
typedef vector<int> Poly;
int F1[maxn],F2[maxn];
int rev[maxn];
void NTT(int * A,int lim,int opt) {
 	int i, j, k, m, gn, g, tmp;
 	for(int i = 0; i < lim; ++i)rev[i] = (i & 1)*(lim >> 1) + (rev[i >> 1] >> 1);
 	for(i = 0; i < lim; ++i)if (rev[i] < i) swap(A[i], A[rev[i]]);
 	for(m = 2; m <= lim; m <<= 1) {
	    k = m >> 1;
	    gn = Pow(3,(p - 1) / m);
	    for(i = 0; i < lim; i += m) {
	    	g = 1;
	    	for (j = 0; j < k; ++j, g = 1ll * g * gn % p) {
	    	    tmp = 1ll * A[i + j + k] * g % p;
	    	    A[i + j + k] = (A[i + j] - tmp + p) % p;
	     	 	A[i + j] = (A[i + j] + tmp) % p;
	    	}
		}
	}
	if(opt == -1){
    	reverse(A+1,A+lim);
    	int inv = Pow(lim,p-2);
    	for(i = 0; i < lim; ++i) A[i] = 1ll * A[i] * inv % p;
	}
}
Poly mul(Poly A,Poly B){
	int n = A.size(),m = B.size();
	int siz = n + m - 1;
	Poly C(siz);
	if(siz < 64){//小于等于64项直接暴力算
		for(int i = 0;i < n;i++){
			for(int j = 0;j < m;j++)C[i+j] = (C[i+j] + 1LL*A[i]*B[j]%p)%p;
		}
		return C;
	}
	int fsiz = 1;
	while(fsiz <= siz)fsiz *= 2;
	for(int i = 0;i < fsiz;i++)F1[i] = F2[i] = 0;
	for(int i = 0;i < n;i++)F1[i] = A[i];
	for(int i = 0;i < m;i++)F2[i] = B[i];
	NTT(F1,fsiz,1);
	NTT(F2,fsiz,1);
	for(int i = 0;i < fsiz;i++)F1[i] = 1ll*F1[i]*F2[i]%p;
	NTT(F1,fsiz,-1);
	for(int i = 0;i < siz;i++){
		C[i] = F1[i];
	}
	return C;
}

与FFT的区别在点值表示法相乘时, NTT需要取模. 另外注意原根需与模数对应

1
for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*B[i]%p;

分治FFT/NTT

设有\(n\)个多项式, 第\(i\)个为\(Func[i]\)

就嗯分治完事了. 时间复杂度\(O(n\log^2 n)\).

1
2
3
4
5
6
vector<Poly>Func;
Poly cal(int L,int R){
	if(R - 1 == L)return Func[L];
	int M = (L+R)/2;
	return mul(cal(L,M),cal(M,R));
}

求逆

板子来源: Siyuan -「算法笔记」多项式求逆

调用invP(Poly A),返回模\(x^n\)意义下\(A\)的逆. 其中\(A\)是一个\(n-1\)次多项式(实际上是一个长度为\(n\)的vector<int>)

时间复杂度: \(O(nlogn)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int P[maxn],Q[maxn];
Poly invP(Poly A){
	int n = A.size();
	int siz = 1;
	while(siz < n)siz *= 2;
	A.resize(siz);
	Poly I(siz);
	I[0] = Pow(A[0],p-2);//inv(A[0]);
	for(int h = 2;h <= siz;h *= 2){
		int H = h*2;
		for(int i = 0;i <= H;i++)P[i] = Q[i] = 0;//重要,必须清空
		for(int i = 0;i < h;i++)P[i] = A[i],Q[i] = I[i];
		NTT(P,H,1);
		NTT(Q,H,1);
		for(int i = 0;i < H;i++){
			P[i] = 1ll*Q[i]*(2 - 1ll*P[i]*Q[i]%p + p)%p;
		}
		NTT(P,H,-1);
		for(int i = 0;i < h;i++)I[i] = P[i];
	}
	I.resize(n);
	return I;
}

求导, 积分

derP求导, intP积分, 复杂度均为\(O(n)\). 其中积分需要先预处理逆元inv[]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
Poly derP(Poly A,bool samelen = 1){
	int n = A.size();
	Poly D(n-1);
	for(int i = 1;i < n;i++)D[i-1] = 1ll*i*A[i]%p;
	if(samelen)D.resize(n);
	return D;
}
Poly intP(Poly A,bool samelen = 1){
	int n = A.size();
	Poly ans(n+1);
	for(int i = 1;i <= n;i++)ans[i] = 1ll*inv[i]*A[i-1]%p;
	if(samelen)ans.resize(n);
	return ans;
}

取对数

输入\(A(x)\), 求\(B(x) \equiv ln(A(x)) \pmod {x^n}\). 两边同时求导得: \[ B'(x) = \frac{A'(x)}{A(x)}\text dx \] 再积分: \[ B(x) = ln(A(x)) = \int B'(x)\text dx = \int\frac{A'(x)}{A(x)}\text dx \] 需要求逆,求导和积分, 需保证\(A[0] = 1\). 复杂度\(O(n\log n)\)

1
2
3
4
5
Poly lnP(Poly A){
	Poly tans = intP(mul(derP(A),invP(A)));
	tans.resize((int)A.size());
	return tans;
}

指数函数

前置: 求\(\ln\)

和求逆基本相同, 返回\(e^{A(x)}\). 复杂度\(O(nlogn)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
Poly expP(Poly A){
	int n = A.size();
	int siz = 1;
	while(siz < n)siz *= 2;
	A.resize(siz);
	Poly E(siz);
	E[0] = 1;
	for(int h = 2;h <= siz;h *= 2){
		Poly B(h),C(h);;
		for(int i = 0;i < h;i++)B[i] = C[i] = E[i];
		B = lnP(B);
		for(int i = 0;i < h;i++)B[i] = (1ll*-B[i] + A[i] + 2ll*p)%p;
		B[0] = (B[0] + 1)%p;
		B = mul(B,C);
		for(int i = 0;i < h;i++)E[i] = B[i];
	}
	E.resize(n);
	return E;
}

快速幂

求\(A(x)^k\), 显然有\(A(x)^k = e^{k\cdot ln(A(x))}\)

保证\(A[0] = 1\)时:

1
2
3
4
5
Poly powP(Poly A,int k){
	A = lnP(A);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	return expP(A);
}

\(A[0] \neq 1\)时:

首先判断是否有\(k \ge n\)且\(A[0] == 0\). 若是, 则返回一个全为\(0\)的多项式. 注意, 必须在取模前判断大小, 且下面的模板没有包含这部分内容

然后对指数取模, 若输入的指数\(k > p\), 则令k1 = k % p, k2 = k % (p - 1) (即\(k2\)对\(\varphi(p)\)取模)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
Poly powP(Poly A,int k1,int k2){
	int n = A.size(),d = 0;
	while(d < n and !A[d])d++;
	if(1ll*d*k1 >= n)return Poly(n);
	A.erase(A.begin(),A.begin()+d);
	int k = Pow(A[0],p-2);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	A = powP(A,k1);
	A.resize(n);
	k = Pow(k,p - 1 - k2);
	for(int i = 0;i < (int)A.size();i++)A[i] = 1ll*A[i]*k%p;
	d *= k1;
	for(int i = n-1;i >= 0;i--){
		A[i] = (i >= d?A[i-d]:0);
	}	
	return A;
}

高阶差分/前缀和

求\(A[]\)的\(k\)阶差分/前缀和. 设\(A\)的OGF: \(F(x) = \sum_{i = 0}^{\infin}A_ix^i\)

显然前缀和是\(F(x)\cdot\sum_{i=0}^{\infin}x^i = F(x)\cdot \frac{1}{1-x}\), 差分是\(F(x)\cdot(1-x)\)

那么\(k\)阶前缀和是\(F(x)\cdot (\frac{1}{1-x})^k\), 差分是\(F(x)\cdot(1-x)^k\). 复杂度为\(O(nlogk)\), 但常数巨大且很不好写.

另有结论:

前缀和\(F(x)\cdot\sum_{i = 0}^nC_{k+i-1}^ix^i\)

差分\(F(x)\cdot\sum_{i=0}^n(-1)^iC_n^ix^i\)

复杂度\(O(n\log n)\). 使用前先预处理逆元.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
Poly presumP(Poly A,ll k){
	int n = A.size();
	Poly G(n);
	G[0] = 1;
	for(int i = 1;i < n;i++)G[i] = G[i-1]*(k+i-1)%p*inv[i]%p;
	A = mul(A,G);
	return A;
}
Poly diffP(Poly A,ll k){
	int n = A.size();
	Poly G(n);
	G[0] = 1;
	for(int i = 1;i < n;i++)G[i] = (G[i-1]*(-1) + p)%p*(k-i+1)%p*inv[i]%p;	
	A = mul(A,G);
	return A;
}

前缀和与差分

区间加多项式

给出\(f(x)\)和\([L,R]\), 对\(a[L...R]\)分别加上\(f(1),f(2),...,f(R-L+1)\).

结论: \(k\)阶多项式的\(k+1\)阶差分余项为常数. 对\(k+1\)阶差分后的余项做\(k+1\)次前缀和即可还原原多项式.

为了维护该操作, 我们先对原数组做\(k+1\)阶差分, 然后算出\(f(1),...,f(k+1)\)的\(k+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
45
46
47
48
49
50
51
52
int d[maxn];
void diff(vector<int>&a,int cnt){
	while(cnt--)for(int i = (int)a.size() - 1;i > 0;i--){a[i] -=a[i-1];if(a[i] < 0)a[i] += p;}
}
void pre(vector<int>&a,int cnt){
	while(cnt--)for(int i = 1;i < (int)a.size();i++)a[i] = (a[i] + a[i-1])%p;
}
ll f(ll x,vector<ll>fac){
	ll res = 0,times = 1;
	for(ll v:fac){
		res = (res + v*times%p)%p;
		times = times*x%p;
	}
	return res;
}
ll g(ll x,vector<ll>fac,ll len){
	return (-f(x+len,fac) + p)%p;
}
signed main(){
	ios::sync_with_stdio(0);
	cin.tie(0);
	cout.setf(ios::fixed,ios::floatfield);
	cout.precision(12);	 
	int n,m,q;
	cin >> n >> m >> q;
	vector<int>a(n+100);
	int num = 6;
	for(int i = 0;i < n;i++)cin >> a[i];
	diff(a,num);
	while(m--){
		int l,r,k;
		cin >> l >> r >> k;
		l--,r--;
		vector<ll>fac(k+1);
		for(int i = k;i >= 0;i--)cin >> fac[i];
		vector<int>add(10+1),sub(10+1);
		for(int i = 0;i <= 10;i++)add[i] = f(i+1,fac),sub[i] = g(i+1,fac,r - l + 1);
		diff(add,num);
		diff(sub,num);
		for(int i = 0;i <= 10;i++)a[i+l] += add[i],a[i+r+1] += sub[i],a[i+l] %= p,a[i+r+1] %= p;
	}
	pre(a,num+1);
	while(q--){
		int l,r;
		cin >> l >> r;
		l--,r--;
		ll ans = (a[r] - (l?a[l-1]:0) + p)%p;
		cout << ans << endl;
	}
	return 0;
}

组合数学

杂乱的结论

长度为n, 可选[1,m]的本质不同序列个数

(本质不同即将序列排序后两个序列不相同) \[ C_{n+m-1}^{m-1} \]

递推组合数

\[ C_n^0 = 1\\ C_n^i = C_n^{i-1}\cdot \frac{n-i+1}{i}\\ C_{n+i-1}^i = C_{n+i-2}^{i-1}\cdot \frac{n+i-1}{i} \]

常用公式

$$ \[\begin{align} \sum_{i = 0}^n C_n^i &= 2^n\\ \sum_{i = 0}^n C_n^i \cdot x^i &= (x+1)^n(二项式定理)\\ \sum_{i = 0,i += 2}^n C_n^i &= 2^{n-1}\\ \sum_{i=1}^nC_n^i \cdot i &= n \cdot 2^{n-1}\\ \sum_{i = 0}^r C_m^i \cdot C_n^{r - i} &= C_{n+m}^r\\ \sum_{i = 0}^n (C_n^i)^2 = C_{2n}^n\\ \end{align}\] $$

常用结论

\[ \sum_{i = 1}^ni*m^{n-i} *C(n,i) = n(m+1)^{n-1} \]

生成函数

OGF

\[ f(x) = \sum_{i = 0}^{\infin}h_ix^i \]

常用公式

  • Part1

\[ \begin{align} \sum_{i=0}^nx^i &= \frac{1-x^{n+1}}{1-x} (等比数列求和)\\ \sum_{i=0}^\infin x^i &= \frac{1}{1-x}\\ \sum_{i = 0}^\infin x^{c\cdot i} &= \frac{1}{1 - x^c}\\ \sum_{i = 0}^{\infin}(i+1)\cdot x^i &= (\frac{1}{1-x})^2 \end{align} \]

  • Part2 \[ \frac{1}{(1-x)^n} = \sum_{i = 0}^\infin C_{n+i-1}^{i-1}x^i\\ \frac{1}{(1-x)^{n+1}} = \sum_{i = 0}^\infin C_{n+i}^i x^i \]

EGF

\[ g(x) = \sum_{i=0}^{\infin}h_x\cdot\frac{x^i}{i!} \]

常用公式

\[ \begin{align} \sum_{i = 0}^{\infin}\frac{x^i}{i!} & = e^x\\ \sum_{i = 0}^{\infin} (-1)^i \frac{x^i}{i!} & = e^{-x} \\ \sum_{i = 0,i +=2}^{\infin}\frac{x^i}{i!} &= \frac{e^x + e^{-x}}{2}\\ \sum_{i = 1,i += 2}^{\infin}\frac{x^i}{i!} & = \frac{e^x - e^{-z}}{2}\\ \sum_{i = 1}^{\infin}(-1)^{i+1} \frac{x^i}{i} & = \ln{(1+x)} \end{align} \]

反演

整除分块

设\(f(x) = \lfloor \frac{k}{x} \rfloor\)

\(f(x)\)分布呈块状,对于任意一个\(i\),有最大的\(j = \lfloor \frac{k}{\lfloor \frac{k}{i} \rfloor} \rfloor\),使得\(f(i) == f(i+1) == ... = f(j)\)

对于类似\(\sum_{i = 1}^{n}\lfloor \frac{k}{i} \rfloor\)的式子,可分块在\(O(\sqrt{n})\)时间内计算完毕

代码:

1
2
3
4
5
for(int l = 1;l <= n;l = r+1){
		if(k/l)r = min(n,k/(k/l));
		else r = n;
		ans += (r-l+1)*(k/l);
	}

对于向上取整, 有 \[ \lceil \frac{k}{x} \rceil = \lfloor \frac{k - 1}{x}\rfloor + 1 \]

狄利克雷(Dirichlet)卷积

定义

设\(f,g\)为两个数论函数,其狄利克雷卷积\(f*g\)为: \[ f \ast g(n) = \sum_{d|n}f(d)g(\frac{n}{d}) \]

常用式子

$$ \[\begin{aligned} \varepsilon=\mu \ast 1&\iff\varepsilon(n)=\sum_{d\mid n}\mu(d)\\ d=1 \ast 1&\iff d(n)=\sum_{d\mid n}1\\ \sigma=\text{id} \ast 1&\iff\sigma(n)=\sum_{d\mid n}d\\ \varphi=\mu \ast \text{id}&\iff\varphi(n)=\sum_{d\mid n}d\cdot\mu(\frac{n}{d}) \\ \end{aligned}\]

$$

常用结论

\[ \sum_{d|n}^n\varphi(d) = n \]

\[ [\gcd(i,j) == 1] = \epsilon(\gcd(i,j) = \sum_{d|\gcd(i,j)}\mu(d) \]

套路题

\(\sum_{i=1}^ngcd(i,n)\)

\[ \begin{align} &\sum_{i=1}^n\gcd(i,n)\\ &=\sum_{d|n}d\sum_{i=1}^n[\gcd(i,n)==d]\\ &=\sum_{d|n}d\sum_{i=1}^{\frac{n}{d}}[\gcd(i,n)==1]\\ &=\sum_{d|n}d\cdotφ(\frac{n}{d}) \end{align} \]

\(\sum_{i = 1}^n\sum_{j=1}^n\gcd(i,j)\)

\[ \begin{align} &\sum_{i = 1}^n\sum_{j=1}^n\gcd(i,j)\\ &= \sum_{d=1}^nd\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j) == d]\\ &=\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[gcd(i,j)==1]\\ &=\sum_{d=1}^nd(2\cdot\sum_{i=1}^{\frac{n}{d}}\varphi(i) - 1) \end{align} \]

关于最后一步,我们设\(f(x) = \sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==1]\)

手画一下,容易看出\(f(x) = 2\cdot\sum_{i=1}^n\varphi(i) - 1\)

另一种做法:

我们知道\(\varphi\)有这么一个性质,\(\sum_{d|n}^{n}\varphi(d) = n\)

也就是\(\sum_{d|\gcd(i,j)}^{n}\varphi(d) = \gcd(i,j)\)

那么: \[ \begin{align} &\sum_{i = 1}^n\sum_{j = 1}^n\gcd(i,j)\\ &=\sum_{i=1}^n\sum_{j=1}^n\sum_{d|\gcd(i,j)}\varphi(d)\\ &=\sum_{d = 1}^n\varphi(d)\sum_{i=1}^n[d|i]\sum_{j=1}^n[d|j]\\ &=\sum_{d = 1}^n\varphi(d)\cdot\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor \end{align} \] 少处理个前缀和=-=而且我感觉这个做法更好理解

例题:洛谷P2398

\(\sum_{i = 1}^{n}\sum_{j=i+1}^{n}\gcd(i,j)\)

和上面那个基本一样 \[ \begin{align} &\sum_{i = 1}^n\sum_{j=i+1}^n\gcd(i,j)\\ &= \sum_{d=1}^nd\sum_{i=1}^n\sum_{j=i+1}^n[gcd(i,j) == d]\\ &=\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\sum_{j=i+1}^{\frac{n}{d}}[gcd(i,j)==1]\\ &=\sum_{d=1}^nd(\sum_{i=1}^{\frac{n}{d}}\varphi(i) - 1) \end{align} \] (最后一步的化简和上一题大同小异,手画一下就有了)

\(\sum_{i=1}^{n}\sum_{i=1}^{m}[\gcd(i,j) == 1]\)

前置知识: \[ \epsilon(n)=\sum_{d|n}\mu(d)= \begin{cases} 1&n=1\\ 0&otherwise \end{cases} \]

开始愉快地推式子: \[ \begin{align} &\sum_{i=1}^{n}\sum_{i=1}^{m}[\gcd(i,j) == 1] \\ &=\sum_{i=1}^{n}\sum_{i=1}^{m}\epsilon(\gcd(i,j))\\ &=\sum_{i=1}^{n}\sum_{i=1}^{m}\sum_{d|gcd(i,j)}\mu(d)\\ &=\sum_{d = 1}^{\min(n,m)}\mu(d)\sum_{i = 1}^n[d|i]\sum_{j = 1}^{m}[d|j]\\ &=\sum_{d=1}^{\min(n,m)}\mu(d)\lfloor \frac{n}{d}\rfloor \lfloor \frac{m}{d} \rfloor \end{align} \] 同理,我们有: \[ \sum_{i = 1}^n\sum_{j = 1}^n\sum_{k = 1}^n [\gcd(i,j,k) == 1] = \sum_{d = 1}^n \mu(d)(\lfloor \frac{n}{d}\rfloor)^3 \]

杂项

小学数学

\[ \sum a_i \cdot b_j = \frac{\sum a_i^2 - \sum b_j^2}{2} \]

因为\((a+b)^2 = a^2 + b^2 + 2ab\)

\[ \sum_{i = 1}^n i = \frac{n \cdot(n+1)}{2}\\ \sum_{i = 1}^n i^2 = \frac{n \cdot (n+1) \cdot(2n + 1)}{6}\\ \sum_{i = 1}^n i^3 = (\frac{n \cdot (n+1)}{2})^2 \]

杂乱的知识点

  • \[ \lfloor \frac{a}{bc} \rfloor = \lfloor \frac{\lfloor \frac{a}{b} \rfloor}{c} \rfloor \]

  • \[ a +b = a\oplus b + 2\cdot a \&b = a|b + a \&b \]

  • 设\(c(n)\)为第\(n\)个Catalan数, 则有\(\sum_{k = 0}^{n-1}c(k)\cdot c(n-1-k) = c(n)\)

  • 区间加平方和(范围\([l,n]\))可通过差分进行. 设差分数组为\(d\), 每次对\([l,n]\)加平方和的操作转换为d[l]++,d[l+1]++, 然后求d的3阶前缀和即可单点查询.

  • 对于一个\(n\)个点的竞赛图(有向完全图), 不经过重复边的最长路最大为\((n为偶数): \frac{n\cdot (n-1)}{2} - \frac{n}{2} + 1\);\((n为奇数):\frac{n\cdot (n-1)}{2}\); 最小为\(n-1\)

图兰定理

对于一个有 \(n\)个点的无向图，若其中不存在三个点的环，则边数不超过 \(n^2/4\)

Pick定理

给定一个顶点均为整点(坐标为整数的点)的简单多边形,其面积\(A\)和内部格点数\(I\),边上格点数\(B\)的关系是: \(A = I+B/2 - 1\)

Cayley公式

一个有\(n\)个节点的完全图有\(n^{n-2}\)种不同的生成树

Polya定理1

-仅旋转,染色数和顶点数均为n.欧拉函数优化(抄自Lskkkno1的题解qwq)
1
2
3
4
5
6
7
8
9
ll polya(int n,int k){
	ll tans = 0;
	for(int i = 1;i*i <= n;i++){
		if(n%i)continue;
		tans += (cal_phi(n/i)*pow_mod(n,i-1))%MOD;
		if(i*i != n)tans += cal_phi(i)*pow_mod(n,n/i-1)%MOD;
	}
	return tans%MOD;
}

\[ C_k^k +C_{k+1}^k +...+C_{n}^k（n > k) \]

某些数列

对于数列\(\{1,2,2,3,3,3,4,4,4,4,...\}\)(数\(i\)出现了\(i\))次, 有通项公式 \[ a_n = \frac{1+2\sqrt{2n}}{2} \]

\[ \sum_{i = 1}^\infin i \cdot \frac{1}{a^i} \]

海伦公式

\(a,b,c\)为三角形边长, \(s = \frac{a+b+c}{2}\)

三角形面积为\(S = \sqrt{s(s-a)(s-b)(s-c)}\)

组合数奇偶性

对于\(C_n^k\), 如果\(n\&k == k\)为奇数, 否则偶数.(按位and)

生成2进制恰好包含k位的下一个数

若\(x\)里有\(k\)个1, 则gen_next(x)为包含\(k\)个1的下一位数

1
2
3
4
5
6
7
int gen_next(int x){
	int b = x & -x;
	int t = x + b;
	int c = t & -t;
	int m = (c/b >> 1) - 1;
	return t | m;
}

求区间[n,m]内不含平方因子的数

筛掉\(p^2\)即可. \(p\)是素数且\(p \in [n,m]\)

等比数列

\(a\)为首项, \(q\)为公比

\[ 求和:S_n = a \cdot \frac{1-q^n}{1-q} (q \ne 1) \]

\[ 求积:P_n = a^n \cdot q^{\frac{n(n-1)}{2}} \]

大素数分解

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
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
namespace prime_fac {

    const int S = 8; // 随机算法判定次数，8~10 就够了

// 龟速乘
    long long mult_mod(long long a, long long b, long long c) {
        a %= c, b %= c;
        long long ret = 0;
        long long tmp = a;
        while (b) {
            if (b & 1) {
                ret += tmp;
                if (ret > c) ret -= c;
            }
            tmp <<= 1;
            if (tmp > c) tmp -= c;
            b >>= 1;
        }
        return ret;
    }

    // 快速幂
    long long qow_mod(long long a, long long n, long long _mod) {
        long long ret = 1;
        long long temp = a % _mod;
        while (n) {
            if (n & 1) ret = mult_mod(ret, temp, _mod);
            temp = mult_mod(temp, temp, _mod);
            n >>= 1;
        }
        return ret;
    }

    // 是合数返回true，不一定是合数返回false
    bool check(long long a, long long n, long long x, long long t) {
        long long ret = qow_mod(a, x, n);
        long long last = ret;
        for (int i = 1; i <= t; i++) {
            ret = mult_mod(ret, ret, n);
            if (ret == 1 && last != 1 && last != n - 1) return true;
            last = ret;
        }
        if (ret != 1) return true;
        return false;
    }

    // 是素数返回true，不是返回false
    mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
    bool Miller_Rabin(long long n) {
        if (n < 2) return false;
        if (n == 2) return true;
        if ((n & 1) == 0) return false;
        long long x = n - 1;
        long long t = 0;
        while ((x & 1) == 0) { x >>= 1; t++; }

        for (int i = 0; i < S; i++) {
            long long a = rng() % (n - 1) + 1;
            if (check(a, n, x, t))
                return false;
        }

        return true;
    }

    long long factor[100];// 存质因数
    int tol; // 质因数的个数，0~tol-1

    long long gcd(long long a, long long b) {
        long long t;
        while (b) {
            t = a;
            a = b;
            b = t % b;
        }
        if (a >= 0) return a;
        return -a;
    }

    long long pollard_rho(long long x, long long c) {
        long long i = 1, k = 2;
        long long x0 = rng() % (x - 1) + 1;
        long long y = x0;
        while (1) {
            i++;
            x0 = (mult_mod(x0, x0, x) + c) % x;
            long long d = gcd(y - x0, x);
            if (d != 1 && d != x) return d;
            if (y == x0) return x;
            if (i == k) { y = x0; k += k; }
        }
    }
    // 对n质因数分解，存入factor，k一般设置为107左右
    void findfac(long long n, int k) {
        if (n == 1) return;
        if (Miller_Rabin(n)) {
            factor[tol++] = n;
            return;
        }
        long long p = n;
        int c = k;
        while (p >= n) p = pollard_rho(p, c--);
        findfac(p, k);
        findfac(n / p, k);
    }
    vector<int> fac(long long n) {
        tol = 0;
        vector<int>ret;
        findfac(n, 107);
        for (int i = 0; i < tol; i++)ret.push_back(factor[i]);
        return ret;
    }
}

子集

子集的和与积

  • 对于集合\(S = \{a_i\}，i = 1...n\), 其全部子集中元素和的和为\(2^{n-1}\cdot \sum_{i = 1}^na_i\)

  • 其全部子集中元素积的和为\(\prod_{i = 1}^n(a_i + 1)\)

  • 从中选出大小为\(k\)的全部子集, 元素积的和可用OGF求解, 即 \[ \sum_{S' \in S \and |S'| = k}\space \prod_{p_i\in S'}p_i = [k] \prod_{i = 1}^n (1+a_i) \] NTT分治即可.

枚举子集

  1. 给定集合\(m\), 枚举它的全部子集:

    1
    for(int S = m;S >= 0;S = (S-1)&m)

  2. 枚举大小为\(n\)的全部集合\(m\), 且对于每个\(m\), 枚举其全部子集\(S\). 时间复杂度\(O(3^n)\):

    1
    2
    3
    4
    for(int m = 0;m < (1<<n);m++){
        for(int S = m;S >= 0; S = (S-1)&m)
    }

子集和/高维前缀和

求解如下问题:

子集和: \[ f[S] = \sum_{T \in S}f[T] \] 超集和: \[ f[T] = \sum_{T\in S}f[S] \]

1
2
3
4
5
6
7
8
for(int i = 0; i < n; i++)
    for(int j = 0; j < (1 << n); j++)
        if(j & (1 << i)){
            f[j] += f[j ^ (1 << i)];//子集和
			//或f[j ^ (1 << i)] += f[j];
        }


时间复杂度\(O(n \cdot 2^n)\)

求双阶乘的最后64位

跑起来很快, 大概是\(O(log^2(n))\)级别 \[ (2n - 1)!! = \prod_{i = 1,i += 2}^{2n-1}i \]

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
67
68
69
70
71
72
73
74
75
#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
lld C[64][64], jc[64];

struct poly
{
    lld x[64];

    poly() { for (int i = 0; i < 64; i++) x[i] = 0; }
    poly(lld c1, lld c0) : poly() { x[1] = c1, x[0] = c0; }
    lld operator()(int k) const { return x[k]; }
    lld &operator()(int k) { return x[k]; }

    poly operator*(const poly &b) const
    {
        poly c;
        for (int i = 0; i < 64; i++)
            for (int j = 0; j <= i; j++)
                c(i) += x[j] * b(i-j);
        return c;
    }

    poly moveto(lld n)
    {
        poly ans;
        for (int i = jc[0] = 1; i < 64; i++)
            jc[i] = jc[i-1] * n;
        for (int i = 0; i < 64; i++)
            for (int j = i; j < 64; j++)
                ans(i) += x[j] * jc[j-i] * C[j][i];
        return ans;
    }
};

map<lld, poly> cache;

poly solve(lld n)
{
    if (cache.count(n)) return cache[n];
    poly half = solve(n / 2);
    half = half * half.moveto(n / 2);
    if (n & 1) half = half * poly(2, n*2-1);
    return cache[n] = half;
}

void init()
{
    poly Px(0, 1);
    cache[0] = Px;

    C[0][0] = C[1][0] = C[1][1] = 1;
    for (int n = 2; n < 64; n++)
    {
        C[n][0] = C[n][n] = 1;
        for (int m = 1; m < n; m++)
            C[n][m] = C[n-1][m-1] + C[n-1][m];
    }
}

int main()
{
    init();
    int T; lld n;
    scanf("%d", &T);

    while (T--)
    {
        scanf("%lld", &n);
        printf("%llu\n", solve(n)(0));
    }

    return 0;
}

\[ \sum_{i = 1}^n \sum_{c \in{cnt} \and c < t[i]}\frac{(n - i)!}{(\prod_{c' \in cnt \and c' \neq c}c'!) \cdot (c-1)!} \]