ACM数学

置顶 stavage

已于 2024-06-30 13:52:51 修改

阅读量1k

点赞数 15

文章标签：算法

于 2023-12-18 16:30:29 首次发布

本文链接：https://blog.csdn.net/m0_73635984/article/details/135064833

版权

基础算法

快速乘&快速幂

LL ksc(LL a, LL b, LL p) {
	LL ans = 0;//快速成是加法
	while (b > 0) {
		if (b & 1) {
			ans = (ans + a) % p;
		}
		b >>= 1;
		a = (a + a) % p;
	}
	return ans;
}
LL ksm(LL a,LL b,LL p)//快速幂
{
	LL ans = 1;//快速幂的乘法
	while (b > 0) {
		if (b & 1) {
			ans = ksc(ans, a, p) % p;//用快速乘代替原来的乘法运算
		}
		b >>= 1;
		a = ksc(a, a, p) % p;
	}
	return ans;
}

欧拉筛

切记主函数初始化init

const int N = 1e6 + 5;
int b[N], prime[N];
int cnt;
void init()
{
    memset(b, 1, sizeof(b));
    b[0] = b[1] = 0;
    for (int i = 2; i < N; i++) {
        if (b[i])
            prime[++cnt] = i;
        for (int j = 1; j <= cnt && prime[j] * i < N; j++)//质数的质数倍，保证不做重复运算
        {
            b[prime[j] * i] = 0;//标记为0，即为非质数
            if (i % prime[j] == 0) break;//解释为i的最小素因子是prime[j]，但我猜测此时i==prime[j]
        }
    }
}

差分

一维还原只需要对差分数组进行前缀和，变换只需要在 l 加 c ，r + 1 减 c

二维差分构造

b[i][j]=a[i][j]+a[i-1][j-1]-a[i-1][j]-a[i][j-1];

还原为二维前缀和

for (int i = 1; i <= n; i++)
    {
        for (int j = 1; j <= m; j++)
        {
            a[i][j]=a[i-1][j]+a[i][j-1]-a[i-1][j-1]+b[i][j];
        }
    }

变换为

b[x1][y1]+=c;
b[x2+1][y1]-=c;
b[x1][y2+1]-=c;
b[x2+1][y2+1]+=c;

快读&快写

在使用__int128的时候使用快读和快写，只需把返回值改为__int128

int qread() {
	int x = 0, y = 1;
	char c = getchar();
	while (c < '0' || c > '9') {
		if (c == '-')
			y = -1;
		c = getchar();
	}
	while (c >= '0' && c <= '9') {
		x = x * 10 + c - '0';
		c = getchar();
	}
	return x * y;
}

/*
//读取操作
for (int i = 1; i <= n; i++) {
	arr[i] = readInt();
	r = max(r, arr[i]);
}
*/
void write(int x) {
     if(x<0) putchar('-'),x=-x;
     if(x>9) write(x/10);
     putchar(x%10+'0');
}

矩阵乘法

多个矩阵相乘，其核心为尽可能用少的时间复杂度，可以考虑先算后面再进行前面的运算顺序

且其一般构造结构体进行存取

ypedef long long LL;
typedef struct {
	LL m[MAX][MAX];
}Mit;

矩阵快速幂/乘

一般都用于求递推公式，将算出来的矩阵设为P，然后需要n阶单位矩阵I进行快速幂

结果一般用S接，然后计算

Mit Mul(Mit a, Mit b)//快速乘
{
	Mit c;
	for (int i = 0; i < MAX; i++)
		for (int j = 0; j < MAX; j++) {
			c.m[i][j] = 0;
			for (int k = 0; k < MAX; k++) {
				a.m[i][k] = (a.m[i][k] % Mod + Mod) % Mod;
				b.m[k][j] = (b.m[k][j] % Mod + Mod) % Mod;
				c.m[i][j] += (a.m[i][k] * b.m[k][j]) % Mod;
			}
			c.m[i][j] = (c.m[i][j] % Mod + Mod) % Mod;
		}
	return c;
}
Mit qpow(Mit P, Mit I, int n)//快速幂
{
	Mit m = P, b = I;
	while (n >= 1) {
		if (n & 1) {
			b = Mul(b, m);
		}
		n = n >> 1;
		m = Mul(m, m);
	}
	return b;
}

高精度

加法

vector<int> x,y,z;
string a, b;
int main()
{
	cin >> a >> b;
	for (int i = a.size() - 1; i >= 0; i--) x.push_back(a[i] - '0');
	for (int i = b.size() - 1; i >= 0; i--) y.push_back(b[i] - '0');
	int t = 0;
	for (int i = 0; i < x.size() || i < y.size(); i++) {
		if (i < x.size()) t += x[i];
		if (i < y.size()) t += y[i];
		z.push_back(t % 10);
		t /= 10;
	}
	if(t) z.push_back(t);//加完之后数位多了的情况
	for (int i = z.size() - 1; i >= 0; i--) {//最高位在最后
		if (i != z.size()) cout << z[i];
		else cout << z[i] << endl;
	}
}

减法

vector<int> x,y,z;
string a, b;
bool pd()//判断x，y两个数据大小
{
	if (x.size() != y.size()) return x.size() > y.size();
	for (int i = x.size() - 1; i >= 0; i--) {
		if (x[i] != y[i]) return x[i] > y[i];
	}
	return 1;
}
int main()
{
	cin >> a >> b;
	for (int i = a.size() - 1; i >= 0; i--) x.push_back(a[i] - '0');
	for (int i = b.size() - 1; i >= 0; i--) y.push_back(b[i] - '0');
	int t = 0;
	int tmp = 0;
	if (pd()) {
		for (int i = 0; i < x.size(); i++) {
			t = x[i] - t;
			if (i < y.size()) t -= y[i];
			z.push_back((t + 10) % 10);
			if (t < 0) t = 1;
			else t = 0;
		}
	}
	else {
		tmp = 1;//标记负号
		for (int i = 0; i < y.size(); i++) {
			t = y[i] - t;
			if (i < x.size()) t -= x[i];
			z.push_back((t + 10) % 10);
			if (t < 0) t = 1;
			else t = 0;
		}
	}
	while (z.size() > 1 && z.back() == 0) z.pop_back();//减法出现0，去除多余的0
	if (tmp) cout << "-";
	for (int i = z.size() - 1; i >= 0; i--) {//也是末尾是最高位
		if (i != z.size()) cout << z[i];
		else cout << z[i] << endl;
	}
	return 0;
}

乘法

一个大数乘一个小数

vector<int> x, y, z;
string a, b;
int c;
int main()
{
	cin >> a >> c;
	for (int i = a.size() - 1; i >= 0; i--) x.push_back(a[i] - '0');
	int t = 0;
	for (int i = 0; i < x.size(); i++) {
		t += x[i] * c;
		z.push_back(t % 10);
		t /= 10;
	}
	while (t) {
		z.push_back(t % 10);
		t /= 10;
	}
	while (z.size() > 1 && z.back() == 0) z.pop_back();//乘0的情况
	for (int i = z.size() - 1; i >= 0; i--) {
		if (i != z.size()) cout << z[i];
		else cout << z[i] << endl;
	}
	return 0;
}

除法

//如果只有除法最好使用正序插入，但由于进行四则运算要与其它运算保持一致所以这里使用逆序
vector<int> x, y, z;
string a, b;
int c;
int main()
{
	cin >> a >> c;
	for (int i = a.size() - 1; i >= 0; i--) x.push_back(a[i] - '0');
	int t = 0;
	for (int i = x.size() - 1; i >= 0; i--) {
		t = t * 10 + x[i];
		z.push_back(t / c);
		t %= c;
	}
	reverse(z.begin(), z.end());
	while (z.size() > 1 && z.back() == 0) z.pop_back();
	for (int i = z.size() - 1; i >= 0; i--) {
		if (i != 0) cout << z[i];
		else cout << z[i] << endl;
	}
	cout << t << endl;//t作余数
	return 0;
}

位运算

计算每一位可以用>>运算将数位后移到个位，进行与1位于运算可算出最后一位是否为1

树状数组 lowbit操作返回最后一位1，eg：10->1010，return 10->2

//位运算>>
while (a[i]) {
            ans[i] += a[i] & 1;
            a[i] = a[i] >> 1;
        }
//lowbit
int lowbit(int x)
{
	return x & -x;
}
while (x) {
			x = x - lowbit(x);
			ans++;
		}

离散化


    vector<int> alls; // 存储所有待离散化的值
	// 二分求出x对应的离散化的值
	int find(int x)
	{
		int l = 0, r = alls.size() - 1;
		while (l < r)
		{
			int mid = l + r >> 1;
			if (alls[mid] >= x) r = mid;
			else l = mid + 1;
		}
		return r + 1;
	}
//main主函数
sort(alls.begin(), alls.end()); // 将所有值排序
alls.erase(unique(alls.begin(), alls.end()), alls.end());	// 去掉重复元素

更加易写但难以理解的离散化：

void work(int c[])
{
	for (int i = 1; i <= n; i++) p[i] = i;
	sort(p + 1, p + n + 1, [&](int x, int y) {
		return a[x] < a[y];
		});//对下标进行排序
	/*相当于bool cmp(int x, int y) {
				return a[x] < a[y];
			}*/
	for (int i = 1; i <= n; i++) a[p[i]] = i; //离散化处理
}

二分

整数

模板一：(用于check检验后答案为*******oxxxxxxx（其中o为答案，x为剩余区间）)

bool check(int k)
{
	
}
int main()
{
	int l = 0, r = 1e9;
	while (l < r) {
		int mid = l + r >> 1;
		if (check(mid)) {
			r = mid;
		}
		else {
			l = mid + 1;
		}
	}
	cout << l << endl;//有可能在l视情况定
	return 0;
}

模板二：(用于check检验后答案为xxxxxo*****（其中o为答案，x为剩余区间）)

bool check(int k)
{
	
}
int main()
{
	int l = 0, r = 1e9;
	while (l < r) {
		int mid = l + r + 1 >> 1;
		if (check(mid)) {
			l = mid;
		}
		else {
			r = mid - 1;
		}
	}
	cout << l << endl;//有可能在l视情况定
	return 0;
}

浮点数

double x;
bool check(double k)
{
    return (k * k * k) < x;//k即为二分结果，要求三次方根，即为k的三次方
}
int main()
{
    l = -10000;
    r = 10000;
    while (r - l > 1e-8)//这里1e-8为精度，一般比要求多2
    {
        double mid = (l + r) / 2;
        if (check(mid)) l = mid;
        else r = mid;
    }
    printf("%.6lf\n", l);
}

三分

LL check(int k)
{
	
}
int main()
{
	int l = 0, r = 1e9;
	while (l <= r) {
		int midl = l + (r - l) / 3;
		int midr = r - (r - l) / 3;
		if (check(midl) <= check(midr)) {
			r = midr - 1;
		}
		else l = midl + 1;
	}                            
	cout << (LL)min(check(l), check(r)) << endl;
	return 0;
}

双指针

for (int i = 0, j = 0; i < n; i ++ )
	{
		while (j < i && check(i, j)) j ++ ;
		//check表示题目的判断条件，没必要一定写函数，一般都是一个for里面套一个while，for里面接具体逻辑
	}

归并排序求逆序对数量

int merge_sort(int l, int r)
{
	if (l == r) return 0;
	int ans = 0;
	int mid = l + r >> 1;
	ans = (merge_sort(l, mid) + merge_sort(mid + 1, r)) % MOD;
	int k = 0, i = l, j = mid + 1;
	while (i <= mid && j <= r) {
		if (b[i] <= b[j]) p[++k] = b[i++];
		else {
			ans = (ans + mid - i + 1) % MOD;
			p[++k] = b[j++];
		}
	}
	while (i <= mid) p[++k] = b[i++];
	while (j <= r) p[++k] = b[j++];
	for (int i = l, j = 1; i <= r; i++, j++) b[i] = p[j];
	return ans;
}

数论

区间筛

切记区间筛的范围为[a,b)

typedef long long LL;
const int N = 1e6;
bool b[N + 5];//a到b的素数判断
bool bprime[N + 5];//0到sqrt(b)的素数判断数组
LL prime[N + 5];//a到b的素数数组
LL isprime[N + 5];//0到sqrt(b)的素数数组
LL p;
LL cnt;
void init(LL l, LL r)
{
	memset(b, 1, sizeof(b));
	memset(bprime, 1, sizeof(bprime));
	for (LL i = 2; i * i < r; i++) {
		if (bprime[i]) {
			isprime[++cnt] = i;
			for (LL j = max(2LL, (l + i - 1) / i) * i; j <= r; j += i) {
				b[j - l + 1] = 0;
			//2LL是2的长整型形式
			// max用于把((a+i-1)/i)*i)和最小素数2进行比较，一般是((a+i-1)/i)*i)，但是但a是1就是2了，它保证删除的都是素数的倍数
			// ((a+i-1)/i)*i)当a为1时，如果没有最大值判定，那么j就是1，会把所有数字删掉
			// 右移回到a，b区间进行判断，判断成功左移回来赋值false
			}
		}
		for (int j = 1; j <= cnt && i * isprime[j] <= N; j++) {//正常欧拉筛写法
			bprime[i * isprime[j]] = 0;
			if (i % isprime[j] == 0) break;
		}
	}
	for (LL i = 1; i <= r - l; i++) {//找到了所有a到b的素数，把素数放进素数表中
		if (b[i]) prime[++p] = i + l - 1;//注意这里如果l是1的话会把1算作素数
	}
}

区间计数

一般使用差分数组计算

区间合并

（将带使用区间存入segs数组中）

void merge(vector<PII> &segs)
	{
		vector<PII> res;

		sort(segs.begin(), segs.end());//sort排序pair优先排序左端

		int st = -2e9, ed = -2e9;//边界
		for (auto seg : segs)
			if (ed < seg.first)
			{
				if (st != -2e9) res.push_back({st, ed});
				st = seg.first, ed = seg.second;
			}
			else ed = max(ed, seg.second);

		if (st != -2e9) res.push_back({st, ed});

		segs = res;
	}
/*
typedef pair<int, int> PII;
vector<PII> vis;
int main()
{
    vis.push_back({ l,r });
    merge(vis);
}
*/

约数

将一个数质因子分解为pow(p1,a1)*pow(p2,a2)*...*pow(pk,ak)

约数个数

约数的个数为(a1+1)*(a2+1)*..*(ak+1)

记得头文件，bits中没有

#include <unordered_map>
typedef long long LL;
unordered_map<LL, LL>sva;//内部自带哈希，但其头文件不在万能头文件中
const int MOD = 1e9 + 7;
int main()
{
	int n;
	cin >> n;
	while (n--) {
		int a;
		cin >> a;
		for (int i = 2; i <= a / i; i++) {//质因子分解过程
			while (a % i == 0) {
				a /= i;
				sva[i]++;
			}
		}
		if (a > 1) sva[a]++;
	}
	LL res = 1;
	for (auto tmp : sva) {
		res = (res * (tmp.second + 1)) % MOD;
	}
	cout << res << endl;
	return 0;
}

素数检验

LL Rand() {//决定了程序的性能
    static LL x = (srand((int)time(0)), rand());
    x += 1000003;
    if (x > 1000000007)
        x -= 1000000007;

    return x;
}
bool Witness(LL a, LL n) {
    LL t = 0, u = n - 1;
    while (!(u & 1))
        u >>= 1, t++;
    LL x = ksm(a, u, n), y;

    while (t--) {
        y = x * x % n;
        if (y == 1 && x != 1 && x != n - 1)
            return true;
        x = y;
    }
    return x != 1;
}
bool MillerRabin(LL n, LL s) {//s为检验次数
    if (n == 2 || n == 3 || n == 5)
        return 1;

    if (n % 2 == 0 || n % 3 == 0 || n % 5 == 0 || n == 1)
        return 0;

    while (s--) {
        if (Witness(Rand() % (n - 1) + 1, n))
            return false;
    }
    return true;
}

大部分时候都会是一个大数的检验，这个时候使用python是更好的

import random
import time

def qread():
    x = 0
    y = 1
    c = input()
    while not c.isdigit():
        if c == '-':
            y = -1
        c = input()
    x = int(c) * y
    return x

def write(x):
    if x < 0:
        print('-', end='')
        x = -x
    if x > 9:
        write(x // 10)
    print(x % 10, end='')

def ksc(a, b, c):
    ans = 0
    while b:
        if b & 1:
            ans = (ans + a) % c
        b >>= 1
        a = (a + a) % c
    return ans

def ksm(a, b, c):
    ans = 1
    while b:
        if b & 1:
            ans = ksc(ans, a, c)
        b >>= 1
        a = ksc(a, a, c)
    return ans

def Rand():
    global x
    x += 1000003
    if x > 1000000007:
        x -= 1000000007
    return x

def Witness(a, n):
    t = 0
    u = n - 1
    while not u & 1:
        u >>= 1
        t += 1
    x = ksm(a, u, n)
    while t:
        y = x * x % n
        if y == 1 and x != 1 and x != n - 1:
            return True
        x = y
        t -= 1
    return x != 1

def MillerRabin(n, s):
    if n == 2 or n == 3 or n == 5:
        return True
    if n % 2 == 0 or n % 3 == 0 or n % 5 == 0 or n == 1:
        return False
    while s:
        if Witness(random.randint(1, n - 1), n):
            return False
        s -= 1
    return True

if __name__ == "__main__":
    y=10**30+50
    print(y)
    x = qread()
    write(x)
    print()
    while not MillerRabin(x, 10):
        x += 1
    write(x)

整除分块

求 $\sum_{i=1}^{n} \left \lfloor \frac{n}{i} \right \rfloor$

int  ans = 0;
for (int l = 1, r; l <= n; l = r + 1) {
	r = n / (n / l);
	ans += n / l * (r - l + 1);
}

求 $\sum_{i=1}^{n} f(i)\left \lfloor \frac{n}{i} \right \rfloor$

需要预处理出f(x)的前缀和数组，再进行分块求和

for (int l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans += (sum[r] - sum[l - 1]) * (n / l);
	}

莫比乌斯函数计算

求单个值

int n;
cin >> n;
int ans = 1;
for (int i = 2; i * i <= n; i++) {
    if (n % i == 0) {
        if (n % (i * i) == 0) {
            cout << 0 << endl;
            return 0;
        }
        while (n % i == 0) n /= i;
        ans *= (-1);
    }
}
if (n) ans *= (-1);
cout << ans << endl;

筛法求值

const int N = 1e4 + 5;
const int max1 = 1e4;
int b[N], prime[N], mo[N];
int cnt;
void init()
{
    memset(b, 1, sizeof(b));
    b[0] = b[1] = 0;
    mo[1] = 1;
    for (int i = 2; i <= max1; i++) {
        if (b[i]) {
            prime[++cnt] = i;
            mo[i] = -1;//质数的莫比乌斯函数值为-1
        }
        for (int j = 1; j <= cnt && prime[j] * i <= max1; j++)//质数的质数倍，保证不做重复运算
        {
            b[prime[j] * i] = 0;//标记为0，即为非质数
            if (i % prime[j] == 0) break;//解释为i的最小素因子是prime[j]，但我猜测此时i==prime[j]
            mo[prime[j] * i] = -mo[i];//由于prime[j]是质数且i不是prime[j]整数倍，所以就会多出一个质数，即变为相反
        }
    }
}

和式的变换

和式的变换主要有几个核心公式

1.交换律，结合律

$2. \sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d\mid gcd(i,j))}d=\sum_{i=1}^{n}\sum_{j=1}^{m}\sum_{d=1}^{n}\left [ d \mid i\right] \left [ d\mid j \right ]d$

$3.\sum_{i=1}^{n}\sum_{j=1}^{m}\left [ gcd(i,j)=k \right ]=\sum_{ik=1}^{n}\sum_{jk=1}^{m}\left [ gcd(ik,jk)=k \right ]=\sum_{i=1}^{\frac{n}{k}}\sum_{j=1}^{\frac{m}{k}}\left [ gcd(i,j)=1 \right ]$

$4.\sum_{i=1}^{n}\sum_{j=1}^{m}A(i)B(j)=\sum_{j=1}^{m}\sum_{i=1}^{n}A(i)B(j)$

$5. \sum_{i=1}^{n}\sum_{j=1}^{m}A(i)B(j)=\sum_{i=1}^{n}A(i)\sum_{j=1}^{m}B(j)$

其中[ ]为示性函数，符合条件的为1

常见变换式：

$1.\left [ gcd(i,j)=1 \right ]=\sum_{d\mid gcd(i,j))}\mu (d)$

同时这个式子也是莫比乌斯反演的常用式子

$2.d\mid gcd(i,j)\rightarrow \left [ d\mid i \right ]\left [ d\mid j \right ]$

$3.\sum_{i=1}^{n}\left [ d\mid i \right ]\rightarrow \left \lfloor \frac{n}{d} \right \rfloor$

$4.\sum_{i=1}^{n}\sum_{j=1}^{m}\left [ gcd(i,j)=k \right ]=\sum_{i=1}^{\frac{n}{k}}\sum_{j=1}^{\frac{m}{k}}\left [ gcd(i,j)=1 \right ]$

其中n/k和m/k为下取整但在代码中int会下取整所以没有标明

5. $\sum_{k=1}^{n}\sum_{d=1}^{\frac{n}{k}}\mu (d)\left \lfloor \frac{n}{kd} \right \rfloor\left \lfloor \frac{m}{kd} \right \rfloor\rightarrow \sum_{k=1}^{n}\sum_{T=1}^{n}\mu (\frac{T}{k})\left \lfloor \frac{n}{T} \right \rfloor\left \lfloor \frac{m}{T} \right \rfloor$

大部分不做这个变换而是在向函数传参时传n/k为n

莫比乌斯反演

狄利克雷卷积

$(f*g)(n)=\sum_{d \mid n}^{}f(d)g(\frac{n}{d})=\sum_{d\mid n}^{}f(\frac{n}{d})g(d)$

三个常用函数

1.元函数

$\varepsilon (n)=\left [ n=1 \right ]$

2.常数函数

1(n)=1

3.恒等函数

id(n)=n

常用的卷积关系

$1.\sum_{d\mid n}^{}\mu (d)=\left [ n=1 \right ]\Leftrightarrow \mu *1=\varepsilon$

$2.\sum_{d\mid n}^{}\varphi (d)=n\Leftrightarrow \varphi *1=id$

$3.\sum_{d\mid n}^{}\mu (d)\frac{n}{d}=\varphi (n)\Leftrightarrow \mu * id=\varphi$

莫比乌斯反演

公式为：
$f(n)=\sum_{d \mid n}g(d)\Leftrightarrow g(n)=\sum_{d \mid n}^{}\mu (d)f(\frac{n}{d})$

但是经常使用的为

$\left [ gcd(i,j)=1 \right ]=\sum_{d\mid gcd(i,j))}\mu (d)$

欧拉反演

欧拉反演只是莫比乌斯反演的一种特殊情况，求

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

经过欧拉反演可以得出

$\sum_{i=1}^{n} \sum_{j=1}^{m}\sum_{d=1}^{n}[d\mid i][d\mid j]\varphi (d)$

这比经过莫比乌斯反演得出的式子简洁许多

欧拉反演的式子：
$\sum_{d\mid n}^{n}\varphi (d)=n$

二项式反演

一般在使用二项式反演时，为已知一个的具体表达式，求另一个

设 f(n) 为至少（或至多）选定个的方案数， g(n) 为恰好选定个的方案数

$1.f(n)=\sum_{i=0}^{n}\binom{n}{i}g(i)\Leftrightarrow g(n)=\sum_{i=0}^{n}(-1)^{n-i}\binom{n}{i}f(i)$ 这类是g为恰好，f为至多

$2.f(n)=\sum_{i=0}^{n}(-1)^{i}\binom{n}{i}g(i)\Leftrightarrow g(n)=\sum_{i=0}^{n}(-1)^{i}\binom{n}{i}f(i)$

$3.f(n)=\sum_{i=n}^{m}\binom{i}{n}g(i)\Leftrightarrow g(n)=\sum_{i=n}^{m}(-1)^{i-n}\binom{i}{n}f(i)$ 这类是g为恰好，f为至少

$4.f(n)=\sum_{i=n}^{m}(-1)^{i}\binom{i}{n}g(i)\Leftrightarrow g(n)=\sum_{i=n}^{m}(-1)^{i}\binom{i}{n}f(i)$

$5.f(n,m)=\sum_{i=0}^{n}\sum_{j=0}^{m}\binom{n}{i}\binom{m}{j}g(i,j)\Leftrightarrow g(n,m)=\sum_{i=0}^{n}\sum_{j=0}^{m}(-1)^{n+m-i-j}\binom{n}{i}\binom{m}{j}f(i,j)$

$6.f(n,m)=\sum_{i=0}^{n}\sum_{j=0}^{m}(-1)^{i+j}\binom{n}{i}\binom{m}{j}g(i,j)\Leftrightarrow g(n,m)=\sum_{i=0}^{n}\sum_{j=0}^{m}(-1)^{i+j}\binom{n}{i}\binom{m}{j}f(i,j)$

$7.f(n,m)=\sum_{i=n}^{?}\sum_{j=m}^{?}\binom{i}{n}\binom{j}{m}g(i,j)\Leftrightarrow g(n,m)=\sum_{i=n}^{?}\sum_{j=m}^{?}(-1)^{i+j-n-m}\binom{i}{n}\binom{j}{m}f(i,j)$

$8.f(n,m)=\sum_{i=n}^{?}\sum_{j=m}^{?}(-1)^{i+j}\binom{i}{n}\binom{j}{m}g(i,j)\Leftrightarrow g(n,m)=\sum_{i=n}^{?}\sum_{j=m}^{?}(-1)^{i+j-}\binom{i}{n}\binom{j}{m}f(i,j)$

杜教筛

杜教筛的核心是狄利克雷卷积，运用狄利克雷卷积把积性函数的前缀和等运算优化到小于线性的时间复杂度

主要根据三个常用积性函数，构造g

杜教筛最好在初始化时，初始化2/3大小的前缀和，否则大概T

//特别注意Sphi和Smo均使用map存储，使用int一般使用杜教筛范围都在1e8-1e9间
//欧拉函数前n项和
int getphi(int n)
{
    if (n < N) return sumphi[n];//筛好的前缀和
    if (Sphi.count(n)) return Sphi[n];//记忆化
    int ans = 1LL * n * (n + 1) / 2;//f * g的前缀和
    for (int l = 2, r; l <= n; l = r + 1) {//数论分块
        r = n / (n / l);
        ans -= (r - l + 1) * getphi(n / l);
    }
    return Sphi[n] = ans;//记忆化
}
//莫比乌斯函数前n项和
int getmo(int n)
{
    if (n < N) return summo[n];
    if (Smo.count(n)) return Smo[n];
    int ans = 1LL;
    for (int l = 2, r; l <= n; l = r + 1) {
        r = n / (n / l);
        ans -= (r - l + 1) * getmo(n / l);
    }
    return Smo[n] = ans;
}
//一般形式
int GetSum(int n) { // 算 f 前缀和的函数
    if (n < N) return f_sum(n);//f的前缀和
    if (S.count(n)) return S[n];//记忆化,S为map定义
    int ans = f_g_sum(n); // 算 f * g 的前缀和
    // 以下这个 for 循环是数论分块
    for (int l = 2, r; l <= n; l = r + 1) { // 注意从 2 开始
        r = (n / (n / l));
        ans -= (g_sum(r) - g_sum(l - 1)) * GetSum(n / l);
        // g_sum 是 g 的前缀和
        // 递归 GetSum 求解
    } return S[n] = ans;
}

欧拉函数计算

求单个值

int a;
while (cin >> a) {
	int ans = a;
	for (int i = 2; i * i <= a; i++) {//寻找素因子
		if (a % i == 0) {
			ans = ans / i * (i - 1);//由于使用公式会出现小数，所以对原公式进行变形
			while (a % i == 0) a /= i;
		}
	}
	if (a > 1) ans = ans / a * (a - 1);
	cout << ans << endl;
}

筛法求值

基于欧拉筛

const int N = 1e6 + 5;
int prime[N], b[N], phi[N];
int cnt;
void init()//欧拉筛法
{
	phi[1] = 1;//第一种情况
	memset(b, 1, sizeof(b));
	b[0] = b[1] = 0;
	for (int i = 2; i < N; i++) {
		if (b[i]) {
			prime[++cnt] = i;
			phi[i] = i - 1;//第四种情况
		}
		for (int j = 1; prime[j] * i < N && j <= cnt; j++) {
			b[prime[j] * i] = 0;
			if (i % prime[j] == 0) {
				phi[prime[j] * i] = phi[i] * prime[j];//第二种情况
				break;
			}
			phi[prime[j] * i] = phi[i] * (prime[j] - 1);//第三种情况
		}
	}
}

逆元

逆元常用的有两种计算方法：

扩展欧几里得和费马小定理，其中exgcd在下文中有写。费马小定理前提为a与MOD互质，则逆元为ksm(a,MOD-2,MOD);

GCD和EXGCD

int gcd(int a, int b)
{
	if (b == 0) return a;
	else return(gcd(b, a % b));
}
void exgcd(LL a, LL b, LL& x, LL& y, LL& d)
{
	if (b == 0) {
		x = 1; y = 0; d = a; return;
	}
	exgcd(b, a % b, x, y, d);
	LL t;
	t = x;
	x = y;
	y = t - (a / b) * y;
}
//逆元算法
if (d == 1) cout << (x + b) % b << endl;
else cout << "没有\n";

线性同余方程

一般使用扩展欧几里得算法，对于ax+by=d有解，d当且仅当是gcd(a,b)的倍数

最小正整数解为：

设 t = b / d * x % m

ans = ( t % ( m / d ) + ( m / d ) ) % ( m % d )

//求解ax同余与b（mod m）
int exgcd(int a, int b, int& x, int& y)
{
	if (!b) {
		x = 1, y = 0;
		return a;
	}
	int d = exgcd(b, a % b, y, x);
	y -= a / b * x;
	return d;
}
int main()
{
	int a, b, m;
	int x, y;
	cin >> a >> b >> m;
	int d = exgcd(a, m, x, y);
	if (b % d != 0) cout << "impossible" << endl;
	else {
		cout << b / d * x % m << endl;//x的值，证明没看懂
	}
}

扩展中国剩余定理（线性同余方程组）

int n, M, ans;
int a[N], m[N];
int ksc(int a, int b, int c)
{
    int ans = 0;
    while (b) {
        if (b & 1) ans = (ans + a) % c;
        b >>= 1;
        a = (a + a) % c;
    }
    return ans;
}
int exgcd(int a, int b, int &x, int &y)
{
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}
int excrt()
{
    M = m[1], ans = a[1];
    for (int i = 2; i <= n; i++) {
        int p = ((a[i] - ans) % m[i] + m[i]) % m[i];
        int x = 0, y = 0;
        int d = exgcd(M, m[i], x, y);
        if (p % d) return -1;//无解
        x = ksc(x, p / d, m[i] / d);
        ans += x * M;
        M *= m[i] / d;
        ans = (ans % M + M) % M;
    }
    return (ans % M + M) % M;
}

高斯消元

线性方程组

const int N = 200;
const double eps = 1e-8;
int n;
double a[N][N];
int gauss()
{
	int r, c;//行，列
	for (c = 1, r = 1; c <= n; c++) {
		int t = r;
		for (int i = r; i <= n; i++) {
			if (fabs(a[i][c]) > fabs(a[t][c])) t = i;
		}
		if (fabs(a[t][c]) < eps) continue;//用于下面判断解的情况
		for (int i = c; i <= n + 1; i++) swap(a[t][i], a[r][i]);
		for (int i = n + 1; i >= c; i--) a[r][i] /= a[r][c];
		for (int i = r + 1; i <= n; i++) {
			if (fabs(a[i][c]) > eps) {
				for (int j = n + 1; j >= c; j--) a[i][j] -= a[r][j] * a[i][c];//目的为了把第r的第c个数变为0
			}
		}
		r++;
	}
	if (r <= n) {
		for (int i = r; i <= n; i++) {
			if (fabs(a[i][n + 1]) > eps) return 2;//无解，出现0 = b情况
		}
		return 1;//无穷多解
	}
	for (int i = n; i >= 1; i--) {
		for (int j = i + 1; j <= n; j++) a[i][n + 1] -= a[i][j] * a[j][n + 1];
	}
	return 0;//有单组解
}





/*
for (int i = 1; i <= n; i++) {
		printf("%.2lf\n", a[i][n + 1]);
	}
*/

异或线性方程组

异或形式则需把a数组改为int形式，在处理上改为异或

const int N = 200;
int n;
int a[N][N];
int gauss()
{
	int r, c;//行，列
	for (c = 1, r = 1; c <= n; c++) {
		int t = r;
		for (int i = r; i <= n; i++) {
			if (a[i][c]) t = i;
		}
		if (!a[t][c]) continue;//用于下面判断解的情况
		for (int i = c; i <= n + 1; i++) swap(a[t][i], a[r][i]);
		for (int i = n + 1; i >= c; i--) a[r][i] /= a[r][c];
		for (int i = r + 1; i <= n; i++) {
			if (a[i][c]) {
				for (int j = n + 1; j >= c; j--) a[i][j] ^= a[r][j];//目的为了把第r的第c个数变为0
			}
		}
		r++;
	}
	if (r <= n) {
		for (int i = r; i <= n; i++) {
			if (a[i][n + 1]) return 2;//无解，出现0 = b情况
		}
		return 1;//无穷多解
	}
	for (int i = n; i >= 1; i--) {
		for (int j = i + 1; j <= n; j++) a[i][n + 1] ^= a[i][j] * a[j][n + 1];
	}
	return 0;//有单组解
}



/*
for (int i = 1; i <= n; i++) {
		printf("%d\n", a[i][n + 1]);
	}
*/

多项式

FFT

const double PI = acos(-1);
struct Complex {
    double x, y;    //  实部和虚部 x + yi
    Complex(double _x = 0.0, double _y = 0.0){ x = _x; y = _y; }
    Complex operator - (const Complex& b) const { return Complex(x - b.x, y - b.y); }
    Complex operator + (const Complex& b) const { return Complex(x + b.x, y + b.y); }
    Complex operator * (const Complex& b) const { return Complex(x * b.x - y * b.y, x * b.y + y * b.x); }
};
//  进行FFT和IFFT前的反转变换
void change(Complex y[], int len)
{
    int i, j, k;
    for (i = 1, j = len / 2; i < len - 1; i++) {
        if (i < j) {
            swap(y[i], y[j]);
        }
        k = len / 2;
        while (j >= k) {
            j -= k;
            k /= 2;
        }
        if (j < k) {
            j += k;
        }
    }
    return;
}
//  FFT
// len必须为2^k形式，这也就导致了在主函数中进行fft前需要使用while (len < m + n + 1) len <<= 1;而在输出结果时需要len = m + n + 1;更改为正确的数值
void fft(Complex y[], int len, int on)
{
    change(y, len);
    for (int h = 2; h <= len; h <<= 1) {
        Complex wn(cos(-on * 2 * PI / h), sin(-on * 2 * PI / h));
        for (int j = 0; j < len; j += h) {
            Complex w(1, 0);
            for (int k = j; k < j + h / 2; k++) {
                Complex u = y[k];
                Complex t = w * y[k + h / 2];
                y[k] = u + t;
                y[k + h / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (on == -1) {
        for (int i = 0; i < len; i++) {
            y[i].x /= len;
        }
    }
}

NTT

原根表

//一个n项和一个m项求卷积
typedef long long LL;
const int p = 998244353, G = 3, Gi = 332748118;//这里的Gi是G的除法逆元
const int N = 5000007;
const double PI = acos(-1);

int n, m;
int res, ans[N];
int len = 1;//
int L;//二进制的位数
int RR[N];
LL a[N], b[N];
inline int qread()
{ 
    int x = 0, y = 1;
    char c = getchar();
    while (c < '0' || c > '9') {
        if (c == '-')
            y = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9') {
        x = x * 10 + c - '0';
        c = getchar();
    }
    return x * y;
}
LL ksm(LL a, LL b, LL p)
{
    LL ans = 1;
    while (b) {
        if (b & 1) ans = ans * a % p;
        b >>= 1;
        a = a * a % p;
    }
    return ans;
}
LL inv(LL x) { return ksm(x, p - 2, p); }
void NTT(LL* A, int type)
{
    for (int i = 0; i < len; ++i)
        if (i < RR[i])
            swap(A[i], A[RR[i]]);
    for (int mid = 1; mid < len; mid <<= 1) {
        LL wn = ksm(G, (p - 1) / (mid * 2), p);
        if (type == -1) wn = ksm(wn, p - 2, p);
        //如果超时了上面if这句话删掉，在下面的if(type == -1)里加上下面这个循环
        /*for (int i = 1; i < len / 2; i ++)
        swap(A[i], A[len - i]); */
        for (int j = mid << 1, pos = 0; pos < len; pos += j) {
            LL w = 1;
            for (int k = 0; k < mid; ++k, w = (w * wn) % p) {
                int x = A[pos + k], y = w * A[pos + mid + k] % p;
                A[pos + k] = (x + y) % p;
                A[pos + k + mid] = (x - y + p) % p;

            }
        }
    }
    if (type == -1) {
        LL limit_inv = inv(len);
        for (int i = 0; i < len; ++i)
            A[i] = (A[i] * limit_inv) % p;
    }
}
//多项式乘法
void poly_mul(LL* a, LL* b, int deg)
{
    for (len = 1, L = 0; len <= deg; len <<= 1) L++;
    for (int i = 0; i < len; ++i) {
        RR[i] = (RR[i >> 1] >> 1) | ((i & 1) << (L - 1));
    }
    NTT(a, 1);
    NTT(b, 1);
    for (int i = 0; i < len; ++i) a[i] = a[i] * b[i] % p;
    NTT(a, -1);
}

int main()
{
    n = qread(), m = qread();
    for (int i = 0; i <= n; ++i) a[i] = (qread() + p) % p;//取模
    for (int i = 0; i <= m; ++i) b[i] = (qread() + p) % p;
    poly_mul(a, b, n + m);
    for (int i = 0; i <= n + m; ++i)
        printf("%d ", a[i]);
    return 0;
}

FWT

//给定长度为2^n的两个序列分别求or,and,xor的结果
const int N = 50007, P = 998244353;
int n, m;
 
void add(int& x, int y) {
    (x += y) >= P && (x -= P);
}
void sub(int& x, int y) {
    (x -= y) < 0 && (x += P);
}
struct FWT {
    int extend(int n) {
        int N = 1;
        for (; N < n; N <<= 1);
        return N;
    }
    void FWTor(vector<int>& a, bool rev) {
        int n = a.size();
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                if (!rev) add(a[i + j + m], a[i + j]);
                else sub(a[i + j + m], a[i + j]);
            }
        }
    }
    void FWTand(vector<int>& a, bool rev) {
        int n = a.size();
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                if (!rev) add(a[i + j], a[i + j + m]);
                else sub(a[i + j], a[i + j + m]);
            }
        }
    }
    void FWTxor(vector<int>& a, bool rev) {
        int n = a.size(), inv2 = (P + 1) >> 1;
        for (int l = 2, m = 1; l <= n; l <<= 1, m <<= 1) {
            for (int j = 0; j < n; j += l) for (int i = 0; i < m; i++) {
                int x = a[i + j], y = a[i + j + m];
                if (!rev) {
                    a[i + j] = (x + y) % P;
                    a[i + j + m] = (x - y + P) % P;
                }
                else {
                    a[i + j] = 1LL * (x + y) * inv2 % P;
                    a[i + j + m] = 1LL * (x - y + P) * inv2 % P;
                }
            }
        }
    }
    vector<int> Or(vector<int> a1, vector<int> a2) {
        int n = max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTor(a1, false);
        a2.resize(N), FWTor(a2, false);
        vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTor(A, true);
        return A;
    }
    vector<int> And(vector<int> a1, vector<int> a2) {
        int n = max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTand(a1, false);
        a2.resize(N), FWTand(a2, false);
        vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTand(A, true);
        return A;
    }
    vector<int> Xor(vector<int> a1, vector<int> a2) {
        int n = max(a1.size(), a2.size()), N = extend(n);
        a1.resize(N), FWTxor(a1, false);
        a2.resize(N), FWTxor(a2, false);
        vector<int> A(N);
        for (int i = 0; i < N; i++) A[i] = 1LL * a1[i] * a2[i] % P;
        FWTxor(A, true);
        return A;
    }
} fwt;
 
int main() {
    scanf("%d", &m);
    n = 1 << m;
    vector<int> a1(n), a2(n);
    for (int i = 0; i < n; i++) scanf("%d", &a1[i]);
    for (int i = 0; i < n; i++) scanf("%d", &a2[i]);
    vector<int> A;
    A = fwt.Or(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    A = fwt.And(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    A = fwt.Xor(a1, a2);
    for (int i = 0; i < n; i++) {
        printf("%d%c", A[i], " \n"[i == n - 1]);
    }
    return 0;
}

博弈论

巴什博弈

一堆石子n个，每次最少取一个，最多取m个，先手必胜点为n%(m+1)!=0

NIM博弈

N堆物品，第i堆物品有Ai个。两名玩家轮流行动，每次可以任选一堆，取走任意多个物品，可把一堆取光，但不能不取。取走最后一件物品者获胜，其先手必胜点为所有对取异或不为0

威佐夫博弈

有两堆各若干个物品，两个人轮流从任意一堆中取出至少一个或者同时从两堆中取出同样多的物品，规定每次至少取一个，至多不限，最后取光者胜利。则是奇异局势先手输，

奇异局势局势的第一个值是未在前面出现过的最小的自然数。

继续分析我们会发现，每种奇异局势的第一个值(这里假设第一堆数目小于第二堆的数目)总是等于当前局势的差值乘上1.618，而1.618 = (sqrt(5.0) + 1) / 2

sg函数

const int N = 1e4 + 5;
void getsg()
{
	memset(sg, 0, sizeof(sg));
	for (int i = 1; i <= N; i++) {
		memset(s, 0, sizeof(s));
		for (int j = 2; j <= 24 && dp[j] <= i; j++) {这里的df[i]为每次可以取走的个数
			s[sg[i - dp[j]]] = 1;
		}
		for (int j = 0; j <= N; j++) {
			if (!s[j]) {
				sg[i] = j;
				break;
			}
		}
	}
}

组合数学

组合数计算

数据a，b均小于2000

const int MOD = 1e9 + 7;
const int N = 2e3 + 5;
int c[N][N];
void init()
{
	for (int i = 0; i < N; i++) {
		for (int j = 0; j <= i; j++) {
			if (!j) c[i][j] = 1;
			else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % MOD;
		}
	}
}
int main()
{
	int a,b;
    cin >> a >> b;
	cout << c[a][b] << endl;
	return 0;
}

数据在1e5以内

typedef long long LL;
const int MOD = 1e9 + 7;
const int N = 1e4 + 5;
LL f[N], inf[N];
LL qmi(LL a, LL b, LL m)
{
	LL ans = 1;
	while (b >= 1) {
		if (b & 1)  ans = (ans * a) % m;
		b = b >> 1;
		a = (a * a) % m;
	}
	return ans;
}
int main()
{
	f[0] = inf[0] = 1;
	for (int i = 1; i <= N; i++) {
		f[i] = f[i - 1] * i % MOD;
		inf[i] = inf[i - 1] * qmi(i, MOD - 2, MOD) % MOD;
	}
	int a, b;
	cin >> a >> b;
	cout << f[a] * inf[b] % MOD * inf[a - b] % MOD << endl;
	return 0;
}

数据达到1e18

typedef long long LL;
LL ksm(LL a, LL b, LL c)
{
	LL ans = 1;//最后答案在ans
	while (b > 0) {
		if (b & 1)//b是奇数的话会出现一个a
		{
			ans = ans * a % c;
		}
		b >>= 1;
		a = (a * a) % c;//b变一半所以a就要平方了
	}
	return ans;
}
LL C(LL n, LL m, LL MOD)
{
	if (m > n) return 0;
	m = min(m, n - m);
	LL a = 1, b = 1;
	while (m) {
		a = a * n % MOD;
		b = b * m % MOD;
		n--;
		m--;
	}
	return a * ksm(b, MOD - 2, MOD) % MOD;
}
LL lucas(LL n, LL m, LL p)
{
	if (m == 0) return 1;
	return C(n % p, m % p, p) * lucas(n / p, m / p, p) % p;
}
int main()
{
	LL a, b;
	cin >> a >> b >> m;
	cout << lucas(a, b) << endl;
	return 0;
}

计算几何

自适应辛普森积分

double f(double x)//原函数
{
	return ;//这里写原函数f(x)表达式
}
double simpson(double l, double r) //Simpson公式
{
	double mid = (l + r) / 2;
	return (f(l) + 4 * f(mid) + f(r)) * (r - l) / 6;
}
double asr(double l, double r, double eps, double tmp)
{
	double mid = (l + r) / 2;
	double ll = simpson(l, mid), rr = simpson(mid, r);
	if (fabs(ll + rr - tmp) <= 15 * eps) return ll + rr + (ll + rr - tmp) / 15.0;     //确认精度
	return asr(l, mid, eps / 2, ll) + asr(mid, r, eps / 2, rr);     //精度不够则递归调用
}


/*
    asr(l, r, eps, simpson(l, r))
*/

基础函数及宏定义变量

const int N = 5007, M = 50007, INF = 0x3f3f3f3f;
const double DINF = 1e18, eps = 1e-8;
const double PI = acos(-1);
struct Point {//二维点
    double x, y;
    Point(double x = 0, double y = 0) :x(x), y(y) { }//构造函数
};
typedef Point Vector;
//向量 + 向量 = 向量，点 + 向量 = 向量
Vector operator + (Vector A, Vector B) { return Vector(A.x + B.x, A.y + B.y); }
//点 - 点 = 向量(向量BC = C - B)
Vector operator - (Point A, Point B) { return Vector(A.x - B.x, A.y - B.y); }
//向量 * 数 = 向量
Vector operator * (Vector A, double p) { return Vector(A.x * p, A.y * p); }
//向量 / 数= 向量
Vector operator / (Vector A, double p) { return Vector(A.x / p, A.y / p); }
//点/向量的比较函数
bool operator < (const Point& a, const Point& b) { return a.x < b.x || (a.x == b.x && a.y < b.y); }
struct Line {//直线定义
    Vector v;
    Point p;
    Line() {}
    Line(Vector v, Point p) :v(v), p(p) {}
    Point get_point_in_line(double t) {//返回直线上一点P = p + v * t
        return p + v * t;
    }
};
//求极角//在极坐标系中，平面上任何一点到极点的连线和极轴的夹角叫做极角。
//单位弧度rad
double Polar_angle(Vector A) { return atan2(A.y, A.x); }
//角度转弧度
double torad(double deg) { return deg / 180 * PI; }
//判断相等
int sgn(double x){//也是dcmp
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
//重载等于运算符
bool operator == (const Point& a, const Point& b){return !sgn(a.x - b.x) && !sgn(a.y - b.y);}
//点积(满足交换律)
double Dot(Vector A, Vector B){return A.x * B.x + A.y * B.y;}
//向量的叉积(不满足交换律)
double Cross(Vector A, Vector B){return A.x *  B.y - B.x * A.y;}
//计算两点距离距离
double dis( Vector a,Vector b)
{
    return sqrt((a.x-b.x)*(a.x-b.x)*1.0+(a.y-b.y)*(a.y-b.y));
}
//判断向量bc是不是向ab的逆时针方向（左边）
bool ToLeftTest(Point a, Point b, Point c){
    return Cross(b - a, c - b) > 0;
}
//取模（模长，求长度）
double Length(Vector A){return sqrt(Dot(A, A));}
//计算两向量夹角
double Angle(Vector A, Vector B){return acos(Dot(A, B) / Length(A) / Length(B));}
// 三角形面积叉乘公式
double Area(Point A, Point B, Point C) { return Cross(B - A, C - A) / 2; } 
//计算两向量构成的平行四边形有向面积
double Area2(Point A, Point B, Point C){return Cross(B - A, C - A);}
//计算向量逆时针旋转九十度之后的单位法线（法向量）
Vector Normal(Vector A) {
    double L = Length(A);
    return Vector(-A.y / L, A.x / L);
}
inline Vector Format(const Vector &A) {
	double L = Length(A);
	return Vector(A.x / L,A.y / L);
}
//计算向量逆时针旋转后的向量
Vector Rotate(Vector A, double rad){
    return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));
}
//点绕着 p 点逆时针旋转 angle
Point rotate(Point p, double angle) {
    Point v = (*this) - p;//这里this不能用，为当前点绕p点的当前点
    double c = cos(angle), s = sin(angle);
    return Point(p.x + v.x * c - v.y * s, p.y + v.x * s + v.y * c);
}
/*
Point rotate(Point q, Point p, double angle) {
    Point v = q - p;//这里this不能用，为当前点绕p点的当前点
    double c = cos(angle), s = sin(angle);
    return Point(p.x + v.x * c - v.y * s, p.y + v.x * s + v.y * c);
}
*/
//将点A绕点B顺时针旋转theta(弧度)
inline Point turn_PP(Point a, Point b, double theta){
    double x = (a.x-b.x) * cos(theta) + (a.y-b.y) * sin(theta) + b.x;
    double y = - (a.x-b.x) * sin(theta) + (a.y-b.y) * cos(theta) + b.y;
    return Point(x,y);
}

复数部分

typedef complex<double> Point;
typedef Point Vector;//复数定义向量后，自动拥有构造函数、加减法和数量积
const double eps = 1e-9;
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
double Length(Vector A){
    return abs(A);
}
double Dot(Vector A, Vector B){//conj(a+bi)返回共轭复数a-bi
    return real(conj(A)*B);
}
double Cross(Vector A, Vector B){
    return imag(conj(A)*B);
}
Vector Rotate(Vector A, double rad){
    return A*exp(Point(0, rad));//exp(p)返回以e为底复数的指数
}

二维计算几何基础

点与线

判断点和直线关系relation

1 在左侧
-1 在右侧
0 在直线上 //直线上两点s和e
A, B:直线上一点,C:待判断关系的点

int relation(Point A, Point B, Point C)
{
    // 1 left -1 right 0 in
    int c = sgn(Cross((B - A), (C - A)));
    if(c < 0) return 1;
    else if(c > 0) return -1;
    return 0;
}

计算两直线交点 Get_line_intersection

调用前要确保两直线p + tv 和 Q + tw之间有唯一交点，当且仅当Corss(v, w) != 0;（t是参数）

Point Get_line_intersection(Point P,Vector v,Point Q,Vector w)
{
	Vector u = P - Q;
	double t = Cross(w, u) / Cross(v, w);
	return P + v * t;
}

计算点到直线的距离 Distance_point_to_line

double Distance_point_to_line(Point P, Point A, Point B)
{
    Vector v1 = B - A, v2 = P - A;
    return fabs(Cross(v1, v2) / Length(v1));//如果不取绝对值，那么得到的是有向距离
}

计算点到线段的距离 Distance_point_to_segment

垂线距离或者PA或者PB距离

double Distance_point_to_segment(Point P, Point A, Point B)
{
    if(A == B) return Length(P - A);//（如果重合那么就是两个点之间的距离，直接转成向量求距离即可）
    Vector v1 = B - A, v2 = P - A, v3 = P - B;
    if(sgn(Dot(v1, v2)) < 0) return Length(v2);//A点左边
    if(sgn(Dot(v1, v3)) > 0)return Length(v3);//B点右边
    return fabs(Cross(v1, v2) / Length(v1));//垂线的距离
}

求点在直线上的投影点 Get_line_projection

Point Get_line_projection(Point P, Point A, Point B)
{
    Vector v = B - A;
    return A + v * (Dot(v, P - A) / Dot(v, v));
}

计算点 P 到直线 AB 的垂足 FootPoint

inline Point FootPoint(Point p, Point a, Point b){
    Vector x = p - a, y = p - b, z = b - a;
    double len1 = Dot(x, z) / Length(z), len2 = - 1.0 * Dot(y, z) / Length(z);//分别计算AP,BP在AB,BA上的投影
    return a + z * (len1 / (len1 + len2));//点A加上向量AF
}

计算点到直线的对称点 Symmetry_PL

inline Point Symmetry_PL(Point p, Point a, Point b){
    return p + (FootPoint(p, a, b) - p) * 2; 
    //实际上就是求垂足之后延长一倍得到的向量，与原来的点加起来就行了
}

判断点是否在线段上 OnSegment

bool OnSegment(Point p, Point a1, Point a2){
    return sgn(Cross(a1-p, a2-p)) == 0 && sgn(Dot(a1-p, a2-p)) <= 0;
}

判断两线段是否相交（不算在端点处相交）segment_proper_intersection

bool segment_proper_intersection(Point a1, Point a2, Point b1, Point b2)
{
    double c1 = Cross(a2 - a1, b1 - a1), c2 = Cross(a2 - a1, b2 - a1);
    double c3 = Cross(b2 - b1, a1 - b1), c4 = Cross(b2 - b1, a2 - b1);
    return sgn(c1) * sgn(c2) < 0 && sgn(c3) * sgn(c4) < 0;
}

判断两线段是否相交（包含端点处相交） Segment_proper_intersection

bool Segment_proper_intersection(Point a1, Point a2, Point b1, Point b2){
    double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
    double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
    //if判断控制是否允许线段在端点处相交，根据需要添加
    if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
        bool f1 = OnSegment(b1, a1, a2);
        bool f2 = OnSegment(b2, a1, a2);
        bool f3 = OnSegment(a1, b1, b2);
        bool f4 = OnSegment(a2, b1, b2);
        bool f = (f1|f2|f3|f4);
        return f;
    }
    return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

判断点是否在一条线段上(不含端点) on_segment

bool on_segment(Point P, Point a1, Point a2)
{
    return sgn(Cross(a1-P, a2-P)) == 0 && sgn(Dot(a1 - P, a2 - P)) < 0;
}

两向量的关系 parallel

直线v和直线Line(s, e);

bool parallel(Line v){ 
 	return sgn((e−s)^(v.e−v.s)) == 0; 
}

两直线关系 linecrossline

0 平行
1 重合
2 相交

int linecrossline(Line v){ 
	if((*this).parallel(v)) 
 	return v.relation(s)==3; 
	return 2; 
}

多边形

三角形四心

外心：三边中垂线交点，到三角形三个顶点距离相同（外接圆圆心）
内心：角平分线的交点，到三角形三边的距离相同（内切圆圆心）
垂心：三条垂线的交点
重心：三条中线的交点，到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点

重心

$\left (\frac{ x1+x2+x3}{3} ,\frac{y1+y2+y3}{3}\right )$

垂心 getcircle

inline Circle getcircle(Point A,Point B,Point C){
    Point P1=(A+B)*0.5,P2=(A+C)*0.5;
    Point O=cross_LL(P1,P1+Normal(B-A),P2,P2+Normal(C-A));
    return Circle(O,Len(A-O));
}

凸多边形面积 convex_polygon_area

叉积的几何意义就是三角形有向面积的二倍，所以这里要除以二

double convex_polygon_area(Point* p, int n)
{
    double area = 0;
    for(int i = 1; i <= n - 2; ++ i)
        area += Cross(p[i] - p[0], p[i + 1] - p[0]);
    return area / 2;
    //return fabs(area / 2);//不加的话求的是有向面积，逆时针为负，顺时针为正
}

非凸多边形面积 polyg_on_area

我们叉积求得的三角形面积是有向的，在外面的面积可以正负抵消掉，
因此非凸多边形也适用，可以从任意点出发划分
可以取原点为起点，减少叉乘次数

double polyg_on_area(Point* p, int n)
{
    double area = 0;
    for(int i = 1; i <= n - 2; ++ i)
        area += Cross(p[i] - p[0], p[i + 1] - p[0]);
    return area / 2;
}

判断点在多边形内 is_point_in_polygon

若点在多边形内返回1，在多边形外部返回0，在多边形上返回-1

int is_point_in_polygon(Point p, vector<Point> poly){//待判断的点和该多边形的所有点的合集
    int wn = 0;
    int n = poly.size();
    for(int i = 0; i < n; ++i){
        if(OnSegment(p, poly[i], poly[(i+1)%n])) return -1;
        int k = sgn(Cross(poly[(i+1)%n] - poly[i], p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(poly[(i+1)%n].y - p.y);
        if(k > 0 && d1 <= 0 && d2 > 0) wn++;
        if(k < 0 && d2 <= 0 && d1 > 0) wn--;
    }
    if(wn != 0)
        return 1;
    return 0;
}

判断点在凸多边形内 ToLeftTest

只需要判断点是否在所有边的左边（按逆时针顺序排列的顶点集）可以使用ToLeftTest ，O ( n )

判断折线（向量）bc是不是向（向量）ab的逆时针方向（左边）转向(ToLeftTest)

bool ToLeftTest(Point a, Point b, Point c){
    return Cross(b - a, c - b) > 0;
}

圆

定义

struct Circle
{
	Point c;
	double r;
	Circle(Point c=Point(),double r=0):c(c),r(r){}
	inline Point point(double a)//通过圆心角求坐标
	{ return Point(c.x+cos(a)*r,c.y+sin(a)*r); }
};
//读入圆
inline Circle read_circle()
{
	Circle C;
	scanf("%lf%lf%lf",&C.c.x,&C.c.y,&C.r);
	return C;
}
//xmult函数只在圆中有应用并不常见
double xmult(Point p1,Point p2,Point p0){
	return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}

圆与直线交点 getLineCircleIntersection

int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol) {
    double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
    double e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - C.r * C.r;
    double delta = f * f - 4 * e * g;//判别式
    if (sgn(delta) < 0)//相离
        return 0;
    if (sgn(delta) == 0) {//相切
        t1 = t2 = -f / (2 * e);
        sol.push_back(L.get_point_in_line(t1));//sol存放交点本身
        return 1;
    }
    //相交
    t1 = (-f - sqrt(delta)) / (2 * e);
    sol.push_back(L.get_point_in_line(t1));
    t2 = (-f + sqrt(delta)) / (2 * e);
    sol.push_back(L.get_point_in_line(t2));
    return 2;
}

求两圆交点 get_circle_circle_intersection

两圆相交保存所有交点返回交点个数（至多两个）

int get_circle_circle_intersection(Circle c1, Circle c2, vector<Point>& sol)
{
    double d = Length(c1.c - c2.c);
    if (sgn(d) == 0) {
        if (sgn(c1.r - c2.r) == 0)return -1;//两圆重合
        return 0;
    }
    if (sgn(c1.r + c2.r - d) < 0)return 0;//相离
    if (sgn(fabs(c1.r - c2.r) - d) > 0)return 0;//在另一个圆的内部

    double a = Polar_angle(c2.c - c1.c);//向量c1c2的极角
    double da = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d));
    //c1c2到c1p1的角
    Point p1 = c1.point(a - da), p2 = c1.point(a + da);
    sol.push_back(p1);
    if (p1 == p2)return 1;
    sol.push_back(p2);
    return 2;
}

点到圆的切线 get_tangents

过点p到圆c的切线，v[i]是第i条切线，返回切线的条数

int get_tangents(Point p, Circle C, Vector* v) {
    Vector u = C.c - p;
    double dist = Length(u);
    if (dist < C.r)return 0;//点在内部，没有切线
    else if (sgn(dist - C.r) == 0) {//p在圆上，只有一条切线
        v[0] = Rotate(u, PI / 2);//切线就是垂直嘛
        return 1;
    }
    else {//否则是两条切线
        double ang = asin(C.r / dist);
        v[0] = Rotate(u, -ang);
        v[1] = Rotate(u, +ang);
        return 2;
    }
}

两圆的公切线 get_tangents

两圆的公切线。
根据两圆的圆心距从小到大排列，一共有6种情况。
情况一：两圆完全重合。有无数条公切线。
情况二：两圆内含，没有公共点。没有公切线。
情况三：两圆内切。有1条外公切线。
情况四：两圆相交。有2条外公切线。
情况五：两圆外切。有3条公切线，其中一条内公切线，两条外公切线。
情况六：两圆相离。有4条公切线，其中内公切线两条，外公切线两条。
可以根据圆心距和半径的关系辨别出这6种情况，然后逐一求解。
情况一和情况二没什么需要求的，情况三和情况五中的内公切线都对应
于“过圆上一点求圆的切线”，只需连接圆心和切点，旋转90°后即可知道
切线的方向向量。这样，问题的关键是求出情况四、五中的外公切线和
情况六中的内公切线。

返回切线的条数，-1表示无穷条切线
a[i]和b[i] 分别是第i条切线在圆A和圆B上的切点

int get_tangents(Circle A, Circle B, Point* a, Point* b)
{
    int cnt = 0;
    if(A.r < B.r)swap(A, B), swap(a, b);
    int d2 = (A.x - B.x) * (A.x - B.x) + (A.y - B.y) * (A.y - B.y);
    int rdiff = A.r - B.r;
    int rsum = A.r + B.r;
    if(d2 < rdiff * rdiff)return 0;//内含
    double base = atan2(B.y - A.y, B.x - A.x);
    if(d2 == 0 && A.r == B.r)return -1;//无限多条切线
    if(d2 == rdiff * rdiff){//内切，1条切线
        a[cnt] = A.point(base);
        b[cnt] = B.point(base); cnt ++ ;
        return 1;
    }
    //有外共切线
    double ang = acos((A.r - B.r) / sqrt(d2));
    a[cnt] = A.point(base + ang);b[cnt] = B.point(base + ang);cnt ++;
    a[cnt] = A.point(base - ang);b[cnt] = B.point(base - ang);cnt ++;
    if(d2 == rsum * rsum){//一条内公切线
        a[cnt] = A.point(base);b[cnt] = B.point(PI + base); cnt ++ ;
    }
    else if(d2 > rsum * rsum){//两条内公切线
        double ang = acos((A.r + B.r) / sqrt(d2));
        a[cnt] = A.point(base + ang); b[cnt] = B.point(PI + base + ang); cnt ++ ;
        a[cnt] = A.point(base - ang); b[cnt] = B.point(PI + base - ang); cnt ++ ;
    }
    return cnt;
}

两圆相交面积 AreaOfOverlap

通过计算两个圆相交所构成的两个扇形面积和减去其构成的筝形的面积

double AreaOfOverlap(Point c1, double r1, Point c2, double r2){
    double d = Length(c1 - c2);
    if(r1 + r2 < d + eps)
        return 0.0;
    if(d < fabs(r1 - r2) + eps){
        double r = min(r1, r2);
        return PI*r*r;
    }
    double x = (d*d + r1*r1 - r2*r2)/(2.0*d);
    double p = (r1 + r2 + d)/2.0;
    double t1 = acos(x/r1);
    double t2 = acos((d - x)/r2);
    double s1 = r1*r1*t1;
    double s2 = r2*r2*t2;
    double s3 = 2*sqrt(p*(p - r1)*(p - r2)*(p - d));
    return s1 + s2 - s3;
}

判直线和圆是否相交,包括相切 intersect_line_circle

需要计算点到直线距离函数

int intersect_line_circle(Point c, double r, Point l1, Point l2) {
    return Distance_point_to_line(c, l1, l2) < r + eps;
}

判线段和圆相交,包括端点和相切 intersect_seg_circle

int intersect_seg_circle(Point c, double r, Point l1, Point l2) {
    double t1 = dis(c, l1) - r, t2 = dis(c, l2) - r;
    Point t = c;
    if (t1 < eps || t2 < eps)
        return t1 > -eps || t2 > -eps;
    t.x += l1.y - l2.y;
    t.y += l2.x - l1.x;
    return xmult(l1, c, t) * xmult(l2, c, t) < eps && Distance_point_to_line(c, l1, l2) - r < eps;
}

判圆和圆相交,包括相切 intersect_circle_circle

int intersect_circle_circle(Point c1, double r1, Point c2, double r2) {
    return dis(c1, c2) < r1 + r2 + eps && dis(c1, c2) > fabs(r1 - r2) - eps;
}

计算圆上到点p最近点,如p与圆心重合,返回p本身 dot_to_circle

Point dot_to_circle(Point c, double r, Point p) {
    Point u, v;
    if (dis(p, c) < eps)
        return p;
    u.x = c.x + r * fabs(c.x - p.x) / dis(c, p);
    u.y = c.y + r * fabs(c.y - p.y) / dis(c, p) * ((c.x - p.x) * (c.y - p.y) < 0 ? -1 : 1);
    v.x = c.x - r * fabs(c.x - p.x) / dis(c, p);
    v.y = c.y - r * fabs(c.y - p.y) / dis(c, p) * ((c.x - p.x) * (c.y - p.y) < 0 ? -1 : 1);
    return dis(u, p) < dis(v, p) ? u : v;
}

计算直线与圆的交点,保证直线与圆有交点 intersection_line_circle

传进来的是直线上两点

void intersection_line_circle(Point c, double r, Point l1, Point l2, Point& p1, Point& p2) {
    Point p = c;
    double t;
    p.x += l1.y - l2.y;
    p.y += l2.x - l1.x;
    p = Get_line_intersection(p, c - p, l1, l2 - l1);
    t = sqrt(r * r - dis(p, c) * dis(p, c)) / dis(l1, l2);
    p1.x = p.x + (l2.x - l1.x) * t;
    p1.y = p.y + (l2.y - l1.y) * t;
    p2.x = p.x - (l2.x - l1.x) * t;
    p2.y = p.y - (l2.y - l1.y) * t;
}

计算圆与圆的交点,保证圆与圆有交点,圆心不重合 intersection_circle_circle

void intersection_circle_circle(Point c1, double r1, Point c2, double r2, Point& p1, Point& p2) {
    Point u, v;
    double t;
    t = (1 + (r1 * r1 - r2 * r2) / dis(c1, c2) / dis(c1, c2)) / 2;
    u.x = c1.x + (c2.x - c1.x) * t;
    u.y = c1.y + (c2.y - c1.y) * t;
    v.x = u.x + c1.y - c2.y;
    v.y = u.y - c1.x + c2.x;
    intersection_line_circle(c1, r1, u, v, p1, p2);
}

求圆外一点对圆(o,r)的两个切点 TangentPoint_PC

poi为圆外的那个点

void TangentPoint_PC(Point poi, Point o, double r, Point& result1, Point& result2) {
    double line = sqrt((poi.x - o.x) * (poi.x - o.x) + (poi.y - o.y) * (poi.y - o.y));
    double angle = acos(r / line);
    Point unitvector, lin;
    lin.x = poi.x - o.x;
    lin.y = poi.y - o.y;
    unitvector.x = lin.x / sqrt(lin.x * lin.x + lin.y * lin.y) * r;
    unitvector.y = lin.y / sqrt(lin.x * lin.x + lin.y * lin.y) * r;
    result1 = Rotate(unitvector, -angle);
    result2 = Rotate(unitvector, angle);
    result1.x += o.x;
    result1.y += o.y;
    result2.x += o.x;
    result2.y += o.y;
    return;
}

求外接圆 get_circumcircle

可以解决三点定圆问题

Circle get_circumcircle(Point p1, Point p2, Point p3)
{
    double Bx = p2.x - p1.x, By = p2.y - p1.y;
    double Cx = p3.x - p1.x, Cy = p3.y - p1.y;
    double D = 2 * (Bx * Cy - By * Cx);
    double ansx = (Cy * (Bx * Bx + By * By) - By * (Cx * Cx + Cy * Cy)) / D + p1.x;
    double ansy = (Bx * (Cx * Cx + Cy * Cy) - Cx * (Bx * Bx + By * By)) / D + p1.y;
    Point p(ansx, ansy);
    return Circle(p, Length(p1 - p));
}

求内接圆 get_incircle

Circle get_incircle(Point p1, Point p2, Point p3)
{
    double a = Length(p2 - p3);
    double b = Length(p3 - p1);
    double c = Length(p1 - p2);
    Point p = (p1 * a + p2 * b + p3 * c) / (a + b + c);
    return Circle(p, Distance_point_to_line(p, p1, p2));
}

判断两圆关系 relationcircle

5 相离
4 外切
3 相交
2 内切
1 内含

int relationcircle(Circle v, Circle p) {
    double d = dis(p.c, v.c);
    if (sgn(d - p.r - v.r) > 0)return 5;
    if (sgn(d - p.r - v.r) == 0)return 4;
    double l = fabs(p.r - v.r
    );
    if (sgn(d - p.r - v.r) < 0 && sgn(d - l) > 0)return 3;
    if (sgn(d - l) == 0)return 2;
    if (sgn(d - l) < 0)return 1;
}

求球面上两个点的最短距离

只允许在球面上走，设两个点的距离为d，球半径为r，则距离为2asin(d/2r)r

ans = asin(d / (2 * r)) * r * 2;

最小圆覆盖 min_cover_circle

p为点集，n为数量，注意从0开始

Circle min_cover_circle(Point p[], int n)
{
    random_shuffle(p, p + n);
    Circle c = { p[0],0 };
    for (int i = 1; i < n; i++) {
        if (sgn(dis(p[i], c.c) - c.r) > 0) {
            c = { p[i],0 };
            for (int j = 0; j < i; j++) {
                if (sgn(dis(p[j], c.c) - c.r) > 0) {
                    c.c.x = (p[i].x + p[j].x) / 2;
                    c.c.y = (p[i].y + p[j].y) / 2;
                    c.r = dis(p[j], c.c);
                    for (int k = 0; k < j; k++) {
                        if (sgn(dis(p[k], c.c) - c.r) > 0) c = get_circumcircle(p[i], p[j], p[k]);
                    }
                }
            }
        }
    }
    return c;
}

常用二维计算几何

凸包 ConvexHull

int ConvexHull(Point* p, int n, Point* ch)
{
    sort(p, p + n);
    int m = 0;
    for (int i = 0; i < n; ++i) {//下凸包
        while (m > 1 && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0)m--;
        ch[m++] = p[i];
    }
    int k = m;
    for (int i = n - 2; i >= 0; --i) {//上凸包
        while (m > k && Cross(ch[m - 1] - ch[m - 2], p[i] - ch[m - 2]) <= 0)m--;
        ch[m++] = p[i];
    }
    if (n > 1) m--;
    return m;
}

判断点是否在凸包内

//须保证A[0]=(0,0)
//t是凸包集合中点的数量
int cmp2(const Point &A, const Point &B) { return Cross(A, B) > 0 || (Cross(A, B) == 0 && Dot(A, A) < Dot(B, B)); }
int inHull(Point a, int t, Point *A)
{
    if (Cross(a, A[0]) > 0 || Cross(A[t-1], a) > 0) return 0;
    int ps = lower_bound(A, A + t, a, cmp2) - A - 1;
    return Cross(a - A[ps], A[(ps + 1) % t] - A[ps]) <= 0;
}

//主函数操作
    Point base = p1[0];
    for (int i = 0; i < t; i++) p1[i] = p1[i] - base;
    inHull(p - base, t, p1);//判断点p

旋转卡壳 Rotating_calipers

Point p[N], con[N];
int n;//凸包返回的点数
double Rotating_calipers()
{
    int op = 1;
    double ans = 0;
    for (int i = 0; i < n; ++i) {
        while (H(A[i], A[i + 1], A[p % n + 1]) >= H(A[i], A[i + 1], A[p])) p = (p + 1) % n;
        //H()为所求的接踵点对条件的计算，这里只卡了一个点可以卡多个点，后续根据实际题目写代码
    }
    return ans;
}

求凸包直径

Point p[N], con[N];
int con_num;//凸包返回的点数
double Rotating_calipers()
{
    int op = 1;
    double ans = 0;
    for (int i = 0; i < con_num; ++i) {
        while (Cross((con[i] - con[op]), (con[i + 1] - con[i])) < Cross((con[i] - con[op + 1]), (con[i + 1] - con[i])))
            op = (op + 1) % con_num;
        ans = max(ans, max(dis(con[i], con[op]), dis(con[i + 1], con[op])));
    }
    return ans;
}

最小矩形覆盖

//计算ac在ab上的投影长度
double project(Vector a, Vector b, Vector c)
{
    return Dot(b - a, c - a) / Length(b - a);
}
void ac(int n, Point *p, Point *rec, double &ans)
{
    for (int i = 0, a = 2, b = 1, c = 2; i < n; i++) {
        auto d = p[i], e = p[i + 1];
        while (sgn(Area2(d, e, p[a]) - Area2(d, e, p[a + 1])) < 0) a = (a + 1) % n;
        while (sgn(project(d, e, p[b]) - project(d, e, p[b + 1])) < 0) b = (b + 1) % n;
        if (!i) c = a;
        while (sgn(project(d, e, p[c]) - project(d, e, p[c + 1])) > 0) c = (c + 1) % n;
        auto x = p[a], y = p[b], z = p[c];
        auto h = Area2(d, e, x) / Length(e - d);
        auto w = Dot(y - z, e - d) / Length(e - d);
        if (h * w < ans) {
            ans = h * w;
            rec[0] = d + (e - d) / Length(e - d) * project(d, e, y);
            rec[3] = d + (e - d) / Length(e - d) * project(d, e, z);
            auto u = Normal(e - d);
            rec[1] = rec[0] + u * h;
            rec[2] = rec[3] + u * h;
        }
    }
}
/*
    printf("%.5Lf\n", ans);
    for (int i = 0; i < 4; i++) {
        if (sgn(p[i].x) == 0) p[i].x = 0;
        if (sgn(p[i].y) == 0) p[i].y = 0;
        printf("%.5Lf %.5Lf\n", p[i].x, p[i].y);
    }
*/

闵可夫斯基和

Vector V1[N], V2[N];
inline int Mincowski(Point *P1, int n, Point *P2, int m, Vector *V)
{
    // 【闵可夫斯基和】求两个凸包{P1}, {P2}的向量集合{V} = {P1 + P2}构成的凸包
    for (int i = 0; i < n; ++i) V1[i] = P1[(i + 1) % n] - P1[i];
    for (int i = 0; i < m; ++i) V2[i] = P2[(i + 1) % m] - P2[i];
    int t = 0, i = 0, j = 0;
    V[t++] = P1[0] + P2[0];
    while (i < n && j < m) V[t++] = V[t - 1] + (sgn(Cross(V1[i], V2[j])) > 0 ? V1[i++] : V2[j++]);
    while (i < n) V[t++] = V[t - 1] + V1[i++];
    while (j < m) V[t++] = V[t - 1] + V2[j++];
    return t;
}

半平面交

最好使用scanf和printf进行输入输出

所用直线（有向）

//切记使用方法，半平面交为向量左侧，向Line中赋值时需要使用Line结构体
/*
Line L[100];
L[i]=Line(a,a-b);//前一个为点，后一个为向量，向量方向为箭头坐标减去箭尾坐标
*/
struct Line {
    Point P;
    Vector v;
    double deg;
    Line() {}
    Line(Point P, Vector v) :P(P), v(v) { deg = atan2(v.y, v.x); }
    bool operator < (const Line& L)const {
        return deg < L.deg;
    }
};
//注意这里的Line定义不同于基础函数中Line的定义使用时需删除基础定义

bool on_left(Line L, Point P) {//注意题目是否包含等号
    return sgn(Cross(L.v, P - L.P)) >= 0;
}
Point get_intersection(Line a, Line b)
{
    Vector u = a.P - b.P;
    double t = Cross(b.v, u) / Cross(a.v, b.v);
    return a.P + a.v * t;
}
int half_plane_intersection(Line* L, int n, Point* poly)
{
    sort(L, L + n);
    int first, last;
    Point* p = new Point[n];
    Line* q = new Line[n];
    q[first = last = 0] = L[0];
    for (int i = 1; i < n; ++i) {
        while (first < last && !on_left(L[i], p[last - 1]))last--;
        while (first < last && !on_left(L[i], p[first]))first++;
        q[++last] = L[i];
        if (fabs(Cross(q[last].v, q[last - 1].v)) < eps) {
            last--;
            if (on_left(q[last], L[i].P))q[last] = L[i];
        }
        if (first < last)p[last - 1] = get_intersection(q[last - 1], q[last]);
    }
    while (first < last && !on_left(q[first], p[last - 1]))last--;
    if (last - first <= 1)return 0;
    p[last] = get_intersection(q[last], q[first]);
    int m = 0;
    for (int i = first; i <= last; ++i)poly[m++] = p[i];
    return m;
}

二分半平面交

eg：给出一个凸多边形海岛求，海岛上距离海最远的距离

进行二分半平面交，二分对象为是否有这个点到各个边的距离是mid，二分半平面需要额外使用Normal函数求法向量

struct Line {
    Point P;
    Vector v;
    double deg;
    Line() {}
    Line(Point P, Vector v) :P(P), v(v) { deg = atan2(v.y, v.x); }
    bool operator < (const Line& L)const {
        return deg < L.deg;
    }
};
bool on_left(Line L, Point P) { return Cross(L.v, P - L.P) > 0; }
Point get_intersection(Line a, Line b)
{
    Vector u = a.P - b.P;
    double t = Cross(b.v, u) / Cross(a.v, b.v);
    return a.P + a.v * t;
}
int half_plane_intersection(Line* L, int n, Point* poly)
{
    sort(L, L + n);
    int first, last;
    Point* p = new Point[n];
    Line* q = new Line[n];
    q[first = last = 0] = L[0];
    for (int i = 1; i < n; ++i) {
        while (first < last && !on_left(L[i], p[last - 1]))last--;
        while (first < last && !on_left(L[i], p[first]))first++;
        q[++last] = L[i];
        if (fabs(Cross(q[last].v, q[last - 1].v)) < eps) {
            last--;
            if (on_left(q[last], L[i].P))q[last] = L[i];
        }
        if (first < last)p[last - 1] = get_intersection(q[last - 1], q[last]);
    }
    while (first < last && !on_left(q[first], p[last - 1]))last--;
    if (last - first <= 1)return 0;
    p[last] = get_intersection(q[last], q[first]);
    int m = 0;
    for (int i = first; i <= last; ++i)poly[m++] = p[i];
    return m;
}
Line L[N];
Point poly[N];
Point p[N];
Vector v[N], v1[N];
int main()
{
    int n;
    scanf("%d", &n);
    for (int i = 0; i < n; i++) {
        scanf("%lf %lf", &p[i].x, &p[i].y);
    }
    p[n] = p[0];
    for (int i = 0; i < n; i++) {
        v[i] = p[i + 1] - p[i];
        v1[i] = Normal(v[i]);//单位法向量，即v[i]向左旋转90
    }
    double l = 0, r = 20000;
    while (r - l > eps) {
        double mid = l + (r - l) / 2;
        for (int i = 0; i < n; i++) {
            L[i] = Line(p[i] + v1[i] * mid, v[i]);
        }
        int m = half_plane_intersection(L, n, poly);
        if (!m) r = mid;
        else l = mid;
    }
    printf("%.6lf\n", l);
    return 0;
}

平面最近点对

#include <bits/stdc++.h>
using namespace std;
const int N = 4e5 + 5;
const double DINF = 1e18;
struct Point {
  double x, y;
  Point(double x = 0, double y = 0) : x(x), y(y) {}
};
Point p[N];
int tmp[N], n;
bool cmp(const Point &a, const Point &b)
{
    if (a.x == b.x)return a.y < b.y;
    return a.x < b.x;
}
bool cmps(const int &a, const int &b)
{
    return p[a].y < p[b].y;
}
double dis(Point p, Point q)
{
    return (p.x - q.x) * (p.x - q.x) + (p.y - q.y) * (p.y - q.y);
}
Point q[N];
double merge(int l, int r)
{
    double d = DINF;
    if (l == r) return d;
    int mid = (l + r) / 2;
    double midx = p[mid].x;
    d = min(merge(l, mid), merge(mid + 1, r));
    int p1 = l, p2 = mid + 1, tot = 0;
    while (p1 <= mid || p2 <= r) {
        if (p1 <= mid && (p2 > r || p[p1].y < p[p2].y)) {
            q[++tot] = p[p1++];
        }
        else {
            q[++tot] = p[p2++];
        }
    }
    for (int i = 1; i <= tot; i++) {
        p[l + i - 1] = q[i];
    }
    tot = 0;
    double dd = d;
    d = sqrt(dd);
    for (int i = l; i <= r; i++) {
        if (abs(p[i].x - midx) <= d) q[++tot] = p[i];
    }
    for (int i = 1; i <= tot; i++) {
        for (int j = i - 1; j >= 1 && q[i].y - q[j].y <= d; j--) {
            dd = min(dd, dis(q[i], q[j]));
            d = sqrt(dd);
        }
    }
    return dd;
}
int main()
{
    scanf("%d", &n);
    for (int i = 1; i <= n; ++i)
        scanf("%lf%lf", &p[i].x, &p[i].y);
    sort(p + 1, p + 1 + n, cmp);
    double ans = merge(1, n);
    printf("%.0lf\n", ans);
    return 0;
}

三角剖分 CulArea

double CulArea(Point A, Point B, Circle C)
{
    Vector OA = A - C.c, OB = B - C.c;
    Vector BA = A - B, BC = C.c - B;
    Vector AB = B - A, AC = C.c - A;
    double DOA = Length(OA), DOB = Length(OB), DAB = Length(AB), r = C.r;
    if (sgn(Cross(OA, OB)) == 0) return 0;
    if (sgn(DOA - C.r) < 0 && sgn(DOB - C.r) < 0) return Cross(OA, OB) * 0.5;
    else if (DOB < r && DOA >= r) {
        double x = (Dot(BA, BC) + sqrtl(r * r * DAB * DAB - Cross(BA, BC) * Cross(BA, BC))) / DAB;
        double TS = Cross(OA, OB) * 0.5;
        return asinl(TS * (1 - x / DAB) * 2 / r / DOA) * r * r * 0.5 + TS * x / DAB;
    }
    else if (DOB >= r && DOA < r) {
        double y = (Dot(AB, AC) + sqrtl(r * r * DAB * DAB - Cross(AB, AC) * Cross(AB, AC))) / DAB;
        double TS = Cross(OA, OB) * 0.5;
        return asinl(TS * (1 - y / DAB) * 2 / r / DOB) * r * r * 0.5 + TS * y / DAB;
    }
    else if (fabs(Cross(OA, OB)) >= r * DAB || Dot(AB, AC) <= 0 || Dot(BA, BC) <= 0) {
        if (Dot(OA, OB) < 0) {
            if (Cross(OA, OB) < 0) return (-acosl(-1.0) - asinl(Cross(OA, OB) / DOA / DOB)) * r * r * 0.5;
            else return (acosl(-1.0) - asinl(Cross(OA, OB) / DOA / DOB)) * r * r * 0.5;
        }
        else return asinl(Cross(OA, OB) / DOA / DOB) * r * r * 0.5;
    }
    else {
        double x = (Dot(BA, BC) + sqrtl(r * r * DAB * DAB - Cross(BA, BC) * Cross(BA, BC))) / DAB;
        double y = (Dot(AB, AC) + sqrtl(r * r * DAB * DAB - Cross(AB, AC) * Cross(AB, AC))) / DAB;
        double TS = Cross(OA, OB) * 0.5;
        return (asinl(TS * (1 - x / DAB) * 2 / r / DOA) + asinl(TS * (1 - y / DAB) * 2 / r / DOB)) * r * r * 0.5 + TS * ((x + y) / DAB - 1);
    }
}
/*
    double res = 0;
    for(int i = 0; i < n; i ++){
        res += CulArea(p[i], p[(i + 1) % n], C);
    }
    cout << fabs(res) << endl;
*/

三维计算几何

基础函数

const double EPS = 0.000001;
typedef struct Point3 {
    double x, y, z;
    Point3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
    bool operator == (const Point3& A) const {
        return x==A.x && y==A.y && z==A.z;
    }
};
typedef Point3 Vector3;

struct Line3   //空间直线
{
    Point3 a, b;
};

struct Plane3   //空间平面
{
    Point3 a, b, c;
    Plane3() {}
    Plane3(Point3 a, Point3 b, Point3 c) :
        a(a), b(b), c(c) { }
};

Point3 read_Point3() {
    double x, y, z;
    scanf("%lf%lf%lf", &x, &y, &z);
    return Point3(x, y, z);
}

Vector3 operator + (const Vector3& A, const Vector3& B) {
    return Vector3(A.x + B.x, A.y + B.y, A.z + B.z);
}

Vector3 operator - (const Point3& A, const Point3& B) {
    return Vector3(A.x - B.x, A.y - B.y, A.z - B.z);
}

Vector3 operator * (const Vector3& A, double p) {
    return Vector3(A.x * p, A.y * p, A.z * p);
}

Vector3 operator / (const Vector3& A, double p) {
    return Vector3(A.x / p, A.y / p, A.z / p);
}
double Dot(Vector3 A, Vector3 B) {
    return A.x * B.x + A.y * B.y + A.z * B.z;
}
double Length(Vector3 A) { return sqrt(Dot(A, A)); }
double Angle(Vector3 A, Vector3 B) { return acos(Dot(A, B) / Length(A) / Length(B)); }
//叉积
Vector3 Cross(const Vector3& A, const Vector3& B) {
    return Vector3(A.y * B.z - A.z * B.y, A.z * B.x - A.x * B.z, A.x * B.y - A.y * B.x);
}
// 三角形abc面积的两倍
double Area2(const Point3& A, const Point3& B, const Point3& C) {
    return Length(Cross(B - A, C - A));
}

矢量差 Subt

Point3 Subt( Point3 u, Point3 v )
{
    Point3 ret;
    ret.x = u.x - v.x;
    ret.y = u.y - v.y;
    ret.z = u.z - v.z;
    return ret;
}

平面法向量 NormalVector

Point3 NormalVector( Plane3 s )
{
    return Cross( Subt( s.a, s.b ), Subt( s.b, s.c ) );
}
Point3 NormalVector( Point3 a, Point3 b, Point3 c )
{
    return Cross( Subt( a, b ), Subt( b, c ) );
}

两点距离 TwoPointDistance

double TwoPointDistance( Point3 p1, Point3 p2 )
{
    return sqrt( (p1.x - p2.x)*(p1.x - p2.x) + (p1.y - p2.y)*(p1.y - p2.y) + (p1.z - p2.z)*(p1.z - p2.z) );
}

点p到平面的距离 DistanceToPlane

// 点p到平面p0-n的距离。n必须为单位向量
double DistanceToPlane(const Point3& p, const Point3& p0, const Vector3& n)
{
    return fabs(Dot(p - p0, n)); // 如果不取绝对值，得到的是有向距离
}

点p在平面上的投影 GetPlaneProjection

// 点p在平面p0-n上的投影。n必须为单位向量(如果不是单位向量就/Length(n)嘛)
Point3 GetPlaneProjection(const Point3& p, const Point3& p0, const Vector3& n)
{
    return p - n * Dot(p - p0, n);
}

直线与平面的交点 LinePlaneIntersection

//直线p1-p2 与平面p0-n的交点
Point3 LinePlaneIntersection(Point3 p1, Point3 p2, Point3 p0, Vector3 n)
{
    Vector3 v = p2 - p1;
    double t = (Dot(n, p0 - p1) / Dot(n, p2 - p1)); //分母为0，直线与平面平行或在平面上
    return p1 + v * t; //如果是线段 判断t是否在0~1之间
}

空间直线距离 LineToLine

//空间直线距离,tmp为两直线的公共法向量
double LineToLine(Line3 u, Line3 v, Point3& tmp)
{
    tmp = Cross(Subt(u.a, u.b), Subt(v.a, v.b));
    return fabs(Dot(Subt(u.a, v.a), tmp)) / Length(tmp);
}

点到直线的距离 DistanceToLine

点p到平面AB距离

double DistanceToLine(const Point3& P, const Point3& A, const Point3& B)
{
    Vector3 v1 = B - A, v2 = P - A;
    return Length(Cross(v1, v2)) / Length(v1);
}

点到线段的距离 DistanceToSeg

p到AB

double DistanceToSeg(Point3 p, Point3 a, Point3 b)
{
    if (a == b)
        return Length(p - a);
    Vector3 v1 = b - a, v2 = p - a, v3 = p - b;
    if (Dot(v1, v2) + EPS < 0)
        return Length(v2);
    else {
        if (Dot(v1, v3) - EPS > 0)
            return Length(v3);
        else
            return Length(Cross(v1, v2)) / Length(v1);
    }
}

求异面直线与的公垂线对应的s LineDistance3D

//求异面直线 p1+s*u与p2+t*v的公垂线对应的s 如果平行|重合，返回false,相交true
bool LineDistance3D(Point3 p1, Vector3 u, Point3 p2, Vector3 v, double& s)
{
    double b = Dot(u, u) * Dot(v, v) - Dot(u, v) * Dot(u, v);
    if (abs(b) <= EPS)
        return false;
    double a = Dot(u, v) * Dot(v, p1 - p2) - Dot(v, v) * Dot(u, p1 - p2);
    s = a / b;
    return true;
}

p1和p2是否在线段a-b的同侧 SameSide

bool SameSide(const Point3& p1, const Point3& p2, const Point3& a, const Point3& b) {
    return Dot(Cross(b - a, p1 - a), Cross(b - a, p2 - a)) - EPS >= 0;
}

判断点P是否在三角形中 PointInTri

bool PointInTri(const Point3& P, const Point3& P0, const Point3& P1, const Point3& P2) {
    return SameSide(P, P0, P1, P2) && SameSide(P, P1, P0, P2) && SameSide(P, P2, P0, P1);
}

判断P是否在三角形A, B, C中，并且到三条边的距离都至少为mindist InsideWithMinDistance

bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, double mindist)
{
    if (!PointInTri(P, A, B, C))
        return false;
    if (DistanceToLine(P, A, B) < mindist)
        return false;
    if (DistanceToLine(P, B, C) < mindist)
        return false;
    if (DistanceToLine(P, C, A) < mindist)
        return false;
    return true;
}

判断P是否在凸四边形中，并且到四条边的距离都至少为mindist InsideWithMinDistance

凸四边形ABCD（顺时针或逆时针）

bool InsideWithMinDistance(const Point3& P, const Point3& A, const Point3& B, const Point3& C, const Point3& D, double mindist)
{
    if (!PointInTri(P, A, B, C))
        return false;
    if (!PointInTri(P, C, D, A))
        return false;
    if (DistanceToLine(P, A, B) < mindist)
        return false;
    if (DistanceToLine(P, B, C) < mindist)
        return false;
    if (DistanceToLine(P, C, D) < mindist)
        return false;
    if (DistanceToLine(P, D, A) < mindist)
        return false;
    return true;
}

三角形P0、P1、P2是否和线段AB相交 TriSegIntersection

bool TriSegIntersection(const Point3& P0, const Point3& P1, const Point3& P2, const Point3& A, const Point3& B, Point3& P)
{
    Vector3 n = Cross(P1 - P0, P2 - P0);

    if (abs(Dot(n, B - A)) <= EPS)
        return false;    // 线段A-B和平面P0P1P2平行或共面
    else   // 平面A和直线P1-P2有惟一交点
    {
        double t = Dot(n, P0 - A) / Dot(n, B - A);

        if (t + EPS < 0 || t - 1 - EPS > 0)
            return false;    // 不在线段AB上 
        P = A + (B - A) * t; // 交点
        return PointInTri(P, P0, P1, P2);
    }
}

空间两三角形是否相交 TriTriIntersection

bool TriTriIntersection(Point3* T1, Point3* T2)
{
    Point3 P;

    for (int i = 0; i < 3; i++) {
        if (TriSegIntersection(T1[0], T1[1], T1[2], T2[i], T2[(i + 1) % 3], P))
            return true;
        if (TriSegIntersection(T2[0], T2[1], T2[2], T1[i], T1[(i + 1) % 3], P))
            return true;
    }

    return false;
}

空间两直线上最近点对 SegSegDistance

返回最近距离

double SegSegDistance(Point3 a1, Point3 b1, Point3 a2, Point3 b2, Point3& ans1, Point3& ans2)
{
    Vector3 v1 = (a1 - b1), v2 = (a2 - b2);
    Vector3 N = Cross(v1, v2);
    Vector3 ab = (a1 - a2);
    double ans = Dot(N, ab) / Length(N);
    Point3 p1 = a1, p2 = a2;
    Vector3 d1 = b1 - a1, d2 = b2 - a2;
    double t1, t2;
    t1 = Dot((Cross(p2 - p1, d2)), Cross(d1, d2));
    t2 = Dot((Cross(p2 - p1, d1)), Cross(d1, d2));
    double dd = Length((Cross(d1, d2)));
    t1 /= dd * dd;
    t2 /= dd * dd;
    ans1 = (a1 + (b1 - a1) * t1);
    ans2 = (a2 + (b2 - a2) * t2);
    return fabs(ans);
}

ab，ac，ad的混合积&四面体体积 Volume6

四面体(三角形体)abcd的有向体积的6倍(六面体(正方形体)的体积)

double Volume6(const Point3& A, const Point3& B, const Point3& C, const Point3& D) {
    return Dot(D - A, Cross(B - A, C - A));
}

四面体的重心 Centroid

Point3 Centroid(const Point3& A, const Point3& B, const Point3& C, const Point3& D) {
    return (A + B + C + D) / 4.0;
}

凸多面体的重心

设vi表示第i个四面体的体积，ai表示重心，则多面体的重心为： $\frac{\sum_{i=1}^{n}a_{i}\times v_{i}}{\sum_{i=1}^{n}v_{i}}$

球的体积交and并

double pow2(double x) { return x * x; }
double pow3(double x) { return x * x * x; }
double dis(Point3 x,Point3 y)
{
    return pow2(x.x - y.x) + pow2(x.y - y.y) + pow2(x.z - y.z);
}
double cos(double a, double b, double c) { return (b * b + c * c - a * a) / (2 * b * c); }

double cap(double r, double h) { return PI * (r * 3 - h) * h * h / 3; }

// 2球体积交
double sphere_intersect(Point3 x, double r1, Point3 y, double r2)
{
    double d = dis(x, y);
    // 相离
    if (d >= pow2(r1 + r2))
        return 0;
    // 包含
    if (d <= pow2(r1 - r2))
        return pow3(min(r1, r2)) * 4 * PI / 3;
    // 相交
    double h1 = r1 - r1 * cos(r2, r1, sqrt(d)), h2 = r2 - r2 * cos(r1, r2, sqrt(d));
    return cap(r1, h1) + cap(r2, h2);
}
// 2球体积并
double sphere_union(Point3 x, double r1, Point3 y, double r2)
{
    double d = dis(x, y);
    // 相离
    if (d >= pow2(r1 + r2))
        return (pow3(r1) + pow3(r2)) * 4 * PI / 3;
    // 包含
    if (d <= pow2(r1 - r2))
        return pow3(max(r1, r2)) * 4 * PI / 3;
    // 相交
    double h1 = r1 + r1 * cos(r2, r1, sqrt(d)), h2 = r2 + r2 * cos(r1, r2, sqrt(d));
    return cap(r1, h1) + cap(r2, h2);
}

三维凸包

凸多面体 ConvexPolyhedron

struct ConvexPolyhedron {
    int n;
    vector<Point3> P, P2;
    vector<Face> faces;

    bool read()
    {
        if (scanf("%d", &n) != 1) {
            return false;
        }

        P.resize(n);
        P2.resize(n);

        for (int i = 0; i < n; i++) {
            P[i] = read_Point3();
            P2[i] = add_noise(P[i]);
        }

        faces = CH3D(P2);
        return true;
    }

    //三维凸包重心
    Point3 centroid()
    {
        Point3 C = P[0];
        double totv = 0;
        Point3 tot(0, 0, 0);

        for (int i = 0; i < faces.size(); i++) {
            Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
            double v = -Volume6(p1, p2, p3, C);
            totv += v;
            tot = tot + Centroid(p1, p2, p3, C) * v;
        }

        return tot / totv;
    }
    //凸包重心到表面最近距离
    double mindist(Point3 C)
    {
        double ans = 1e30;

        for (int i = 0; i < faces.size(); i++) {
            Point3 p1 = P[faces[i].v[0]], p2 = P[faces[i].v[1]], p3 = P[faces[i].v[2]];
            ans = min(ans, fabs(-Volume6(p1, p2, p3, C) / Area2(p1, p2, p3)));
        }

        return ans;
    }
};

凸包的定义 Face

struct Face {
    int v[3];
    Face(int a, int b, int c)
    {
        v[0] = a;
        v[1] = b;
        v[2] = c;
    }//逆时针旋转
    Vector3 Normal(const vector<Point3>& P) const
    {
        return Cross(P[v[1]] - P[v[0]], P[v[2]] - P[v[0]]);
    }
    // f是否能看见P[i]
    int CanSee(const vector<Point3>& P, int i) const
    {
        return Dot(P[i] - P[v[0]], Normal(P)) > 0;
    }
};

加干扰防止多点共面 add_noise

double rand01()
{
    return rand() / (double)RAND_MAX;
}
double randeps()
{
    return (rand01() - 0.5) * EPS;
}
Point3 add_noise(const Point3& p)
{
    return Point3(p.x + randeps(), p.y + randeps(), p.z + randeps());
}

增量法求三维凸包 CH3D

// 注意：没有考虑各种特殊情况（如四点共面）。实践中，请在调用前对输入点进行微小扰动
//vector<Face>CH3D(Point_3D* p, int n)//所有面的点集和点数
vector<Face> CH3D(const vector<Point3>& P)
{
    int n = P.size();
    vector<vector<int> > vis(n);
    for (int i = 0; i < n; i++) {
        vis[i].resize(n);
    }
    vector<Face> cur;
    cur.push_back(Face(0, 1, 2)); // 由于已经进行扰动，前三个点不共线
    cur.push_back(Face(2, 1, 0));

    for (int i = 3; i < n; i++) {
        vector<Face> next;
        for (int j = 0; j < cur.size(); j++) {
            Face& f = cur[j];
            int res = f.CanSee(P, i);

            if (!res) {
                next.push_back(f);
            }

            for (int k = 0; k < 3; k++) {
                vis[f.v[k]][f.v[(k + 1) % 3]] = res;
            }
        }

        for (int j = 0; j < cur.size(); j++)
            for (int k = 0; k < 3; k++) {
                int a = cur[j].v[k], b = cur[j].v[(k + 1) % 3];

                if (vis[a][b] != vis[b][a] && vis[a][b]) // (a,b)是分界线，左边对P[i]可见
                {
                    next.push_back(Face(a, b, i));
                }
            }

        cur = next;
    }

    return cur;
}

三维凸包模板合集

const int MAXN=550;
const double eps=1e-8; 
struct Point
{
    double x,y,z;
    Point(){}
 
    Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
 
    //两向量之差
    Point operator -(const Point p1)
    {
        return Point(x-p1.x,y-p1.y,z-p1.z);
    }
 
    //两向量之和
    Point operator +(const Point p1)
    {
        return Point(x+p1.x,y+p1.y,z+p1.z);
    }
 
    //叉乘
    Point operator *(const Point p)
    {
        return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
    }
 
    Point operator *(double d)
    {
        return Point(x*d,y*d,z*d);
    }
 
    Point operator / (double d)
    {
        return Point(x/d,y/d,z/d);
    }
 
    //点乘
    double  operator ^(Point p)
    {
        return (x*p.x+y*p.y+z*p.z);
    }
};
 
struct CH3D
{
    struct face
    {
        //表示凸包一个面上的三个点的编号
        int a,b,c;
        //表示该面是否属于最终凸包上的面
        bool ok;
    };
    //初始顶点数
    int n;
    //初始顶点
    Point P[MAXN];
    //凸包表面的三角形数
    int num;
    //凸包表面的三角形
    face F[8*MAXN];
    //凸包表面的三角形
    int g[MAXN][MAXN];
    //向量长度
    double vlen(Point a)
    {
        return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
    }
    //叉乘
    Point cross(const Point &a,const Point &b,const Point &c)
    {
        return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),
                     (b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),
                     (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)
                     );
    }
    //三角形面积*2
    double area(Point a,Point b,Point c)
    {
        return vlen((b-a)*(c-a));
    }
    //四面体有向体积*6
    double volume(Point a,Point b,Point c,Point d)
    {
        return (b-a)*(c-a)^(d-a);
    }
    //正：点在面同向
    double dblcmp(Point &p,face &f)
    {
        Point m=P[f.b]-P[f.a];
        Point n=P[f.c]-P[f.a];
        Point t=p-P[f.a];
        return (m*n)^t;
    }
    void deal(int p,int a,int b)
    {
        int f=g[a][b];//搜索与该边相邻的另一个平面
        face add;
        if(F[f].ok)
        {
            if(dblcmp(P[p],F[f])>eps)
              dfs(p,f);
            else
            {
                add.a=b;
                add.b=a;
                add.c=p;//这里注意顺序，要成右手系
                add.ok=true;
                g[p][b]=g[a][p]=g[b][a]=num;
                F[num++]=add;
            }
        }
    }
    void dfs(int p,int now)//递归搜索所有应该从凸包内删除的面
    {
         F[now].ok=0;
         deal(p,F[now].b,F[now].a);
         deal(p,F[now].c,F[now].b);
         deal(p,F[now].a,F[now].c);
    }
    bool same(int s,int t)
    {
        Point &a=P[F[s].a];
        Point &b=P[F[s].b];
        Point &c=P[F[s].c];
        return fabs(volume(a,b,c,P[F[t].a]))<eps &&
               fabs(volume(a,b,c,P[F[t].b]))<eps &&
               fabs(volume(a,b,c,P[F[t].c]))<eps;
    }
    //构建三维凸包
    void create()
    {
        int i,j,tmp;
        face add;
 
        num=0;
        if(n<4)return;
    //**********************************************
        //此段是为了保证前四个点不共面
        bool flag=true;
        for(i=1;i<n;i++)
        {
            if(vlen(P[0]-P[i])>eps)
            {
                swap(P[1],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前三个点不共线
        for(i=2;i<n;i++)
        {
            if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps)
            {
                swap(P[2],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
        flag=true;
        //使前四个点不共面
        for(int i=3;i<n;i++)
        {
            if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps)
            {
                swap(P[3],P[i]);
                flag=false;
                break;
            }
        }
        if(flag)return;
    //*****************************************
        for(i=0;i<4;i++)
        {
            add.a=(i+1)%4;
            add.b=(i+2)%4;
            add.c=(i+3)%4;
            add.ok=true;
            if(dblcmp(P[i],add)>0)swap(add.b,add.c);
            g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;
            F[num++]=add;
        }
        for(i=4;i<n;i++)
        {
            for(j=0;j<num;j++)
            {
                if(F[j].ok&&dblcmp(P[i],F[j])>eps)
                {
                    dfs(i,j);
                    break;
                }
            }
        }
        tmp=num;
        for(i=num=0;i<tmp;i++)
          if(F[i].ok)
            F[num++]=F[i];
 
    }
    //表面积
    double area()
    {
        double res=0;
        if(n==3)
        {
            Point p=cross(P[0],P[1],P[2]);
            res=vlen(p)/2.0;
            return res;
        }
        for(int i=0;i<num;i++)
          res+=area(P[F[i].a],P[F[i].b],P[F[i].c]);
        return res/2.0;
    }
    double volume()
    {
        double res=0;
        Point tmp(0,0,0);
        for(int i=0;i<num;i++)
           res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
        return fabs(res/6.0);
    }
    //表面三角形个数
    int triangle()
    {
        return num;
    }
    //表面多边形个数
    int polygon()
    {
        int i,j,res,flag;
        for(i=res=0;i<num;i++)
        {
            flag=1;
            for(j=0;j<i;j++)
              if(same(i,j))
              {
                  flag=0;
                  break;
              }
            res+=flag;
        }
        return res;
    }
    //三维凸包重心
    Point barycenter()
    {
        Point ans(0,0,0),o(0,0,0);
        double all=0;
        for(int i=0;i<num;i++)
        {
            double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
            ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
            all+=vol;
        }
        ans=ans/all;
        return ans;
    }
    //点到面的距离
    double ptoface(Point p,int i)
    {
        return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
    }
};
CH3D hull;
int main()
{
    while(scanf("%d",&hull.n)==1)
    {
        for(int i=0;i<hull.n;i++)
        {
            scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
        }
        hull.create();
        Point p=hull.barycenter();//求重心
        double ans1=1e20;
        for(int i=0;i<hull.num;i++)
        {
            ans1=min(ans1,hull.ptoface(p,i));
        }
        printf("%.3f\n",ans1);
    }
    return 0;
}

参考资料

OI Wiki

繁凡的ACM算法全家桶实在tql