笔记2 (Day6-Day?)

Day6 动态规划优化

单调队列、单调栈优化

P2698 [USACO12MAR] Flowerpot S

二分+滑动窗口 \(\mathcal{O} (n\log n)\)
或者两个单调队列 \(\mathcal{O} (n)\)

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
const int INF=1e9;
int d,n,ans=INF;
struct water{
	int x,y;
	bool operator <(const water& b) const{
		return x==b.x? y<b.y:x<b.x;
	}
}a[N];
deque<int> q1,q2;// q1最小,q2最大
int l,r;

int main(){
	scanf("%d%d",&n,&d);
	for(int i=1;i<=n;i++){
		scanf("%d%d",&a[i].x,&a[i].y);
	}
	sort(a+1,a+n+1);
	q1.push_back(1);q2.push_back(1);
	for(l=1,r=1;l<=n;l++){
		while(!q1.empty()&&q1.front()<l)q1.pop_front();
		while(!q2.empty()&&q2.front()<l)q2.pop_front();
		while(r<n&&a[q2.front()].y-a[q1.front()].y<d){
			r++;
			while(!q1.empty()&&a[q1.back()].y>a[r].y)q1.pop_back();
			while(!q2.empty()&&a[q2.back()].y<a[r].y)q2.pop_back();
			q1.push_back(r);q2.push_back(r);
		}
		if(a[q2.front()].y-a[q1.front()].y>=d){
			ans = min(ans,a[r].x-a[l].x);
		}
	}
	printf("%d\n",ans==INF? -1:ans);
	
	return 0;
}

P3572 [POI 2014] PTA-Little Bird

单调队列优化 \(dp\)

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
int n,d[N],Q,f[N];
deque<int> q;

int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++)scanf("%d",&d[i]);
	scanf("%d",&Q);
	while(Q--){
		int k;scanf("%d",&k);
		memset(f,0,sizeof f);
		if(q.size())q.clear();
		f[1]=0;q.push_back(1);
		for(int i=2;i<=n;i++){
			while(!q.empty()&&i-q.front()>k)q.pop_front();
			f[i]=f[q.front()]+(d[i]>=d[q.front()]);
			while(!q.empty()&&(f[q.back()]>f[i]||(f[q.back()]==f[i]&&d[q.back()]<d[i])))q.pop_back();
			q.push_back(i);
		}
		printf("%d\n",f[n]);
	}
	
	return 0;
}

P3594 [POI 2015 R3] 狼坑 Trous de loup

注意到取长度为 \(d\) 的区间为 \(0\) 一定更优

#include<bits/stdc++.h>
using namespace std;
#define ll long long
const int N=2e6+5;
int n,d,res,l=1,head,tail;
ll p,w[N],sum[N],q[N];

int main(){
	cin>>n>>p>>d;
	for(int i=1;i<=n;++i){
		cin>>w[i];
		sum[i]=sum[i-1]+w[i];
	}
	for(int r=1;r<=n;++r){
		while(head<=tail&&sum[r]-sum[r-d]>=sum[q[tail]]-sum[q[tail]-d])tail--;
		q[++tail]=r;
		while(l<=r&&sum[r]-sum[l-1]-(sum[q[head]]-sum[q[head]-d])>p){
			l++;
			while(head<=tail&&q[head]-d+1<l)head++;
		}
		if(l<=r)res=max(res,r-l+1);
	}
	cout<<res<<'\n';
	return 0;
}

P6510 奶牛排队

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
int n,h[N],ans,q2[N],q1[N],top1,top2;

int main(){
	cin.tie(0);cout.tie(0);ios::sync_with_stdio(0);
	cin>>n;
	for(int i=1;i<=n;i++)
		cin>>h[i];
	for(int i=1;i<=n;i++){
		while(top1&&h[q1[top1]]<h[i])top1--;
		while(top2&&h[q2[top2]]>=h[i])top2--;
		int pos=upper_bound(q2+1,q2+top2+1,q1[top1])-q2;
		if(pos<=top2){
			ans=max(ans,i-q2[pos]+1);
		}
		q2[++top2]=i;q1[++top1]=i;
	}
	cout<<ans<<'\n';
	return 0;
}

斜率优化

我用AI总结的办法:
在斜率优化 DP 中，最让人头疼的莫过于 “凸壳方向” 和 “不等式方向”。网上常见的方法是用几何直觉（上凸壳 / 下凸壳）来判断，但这容易混淆。提供一种 纯代数推导 的方法：通过假设“新决策比旧决策更优”来推出凸壳维护的所有条件，无需记忆任何方向。

---

一、问题标准化

设 DP 转移为：

\[dp_i = \operatorname{opt}\limits_{0 \le j < i} \Big( dp_j + w(j, i) \Big) + C_i \]

其中 \(w(j, i)\) 能展开为关于 \(S_i\)（前缀和）的二次式，通常可整理为：

\[dp_i = \operatorname{opt}\limits_{j < i} \Big( Y_j + K_j \cdot X_i \Big) + D_i \]

这里：

• \(X_i\) 是查询点，单调递增（大部分题目中为前缀和）；
• \(K_j\) 是决策 \(j\) 的斜率，单调递增（若递减，可对 \(X\) 取负或后续调整）；
• \(Y_j\) 是决策 \(j\) 的截距；
• \(\operatorname{opt}\) 为 \(\min\) 或 \(\max\)。

---

二、核心推导：两个决策的比较

考虑两个决策 \(j < k\)（编号顺序）。在某个查询点 \(x = X_i\) 下，若决策 \(k\) 比 \(j\) 更优，则：

• 若求 最大值（\(\operatorname{opt} = \max\)）：

\[Y_k + K_k \cdot x \ge Y_j + K_j \cdot x \]

• 若求 最小值（\(\operatorname{opt} = \min\)）：

\[Y_k + K_k \cdot x \le Y_j + K_j \cdot x \]

两种情形统一可写成：

\[\left( Y_k - Y_j \right) \ge \left( K_j - K_k \right) \cdot x \quad (\text{当 } \operatorname{opt}=\max) \]

\[\left( Y_k - Y_j \right) \le \left( K_j - K_k \right) \cdot x \quad (\text{当 } \operatorname{opt}=\min) \]

由于 \(K\) 单调递增，故 \(K_j - K_k < 0\)。为了方便，定义：

\[\text{slope}(j, k) = \frac{Y_k - Y_j}{K_j - K_k} = - \frac{Y_k - Y_j}{K_k - K_j} \]

但更常用的定义是两点间斜率 \(\frac{Y_k - Y_j}{K_k - K_j}\)，我们稍作变换。

---

2.1 情形一：\(K\) 单调递增，\(\operatorname{opt} = \max\)

由 \(Y_k - Y_j \ge (K_j - K_k) \cdot x\)，由于 \(K_j - K_k < 0\)，两边同除以负数，不等号反向：

\[\frac{Y_k - Y_j}{K_j - K_k} \le x \]

即：

\[x \ge -\frac{Y_k - Y_j}{K_k - K_j} \]

记：

\[\operatorname{limit}(j, k) = -\frac{Y_k - Y_j}{K_k - K_j} \]

则当 \(x \ge \operatorname{limit}(j, k)\) 时，\(k\) 优于 \(j\)；当 \(x < \operatorname{limit}\) 时，\(j\) 优于 \(k\)。

因为 \(x\) 单调递增，一旦 \(x\) 达到这个阈值，\(k\) 就永远优于 \(j\)。所以最优决策会随着 \(x\) 的增大而 向斜率更大的决策移动。

---

2.2 情形二：\(K\) 单调递增，\(\operatorname{opt} = \min\)

由 \(Y_k - Y_j \le (K_j - K_k) \cdot x\)，同样除以负数：

\[\frac{Y_k - Y_j}{K_j - K_k} \ge x \]

即：

\[x \le -\frac{Y_k - Y_j}{K_k - K_j} \]

此时当 \(x \le \operatorname{limit}\) 时，\(k\) 优于 \(j\)，当 \(x\) 超过阈值后，\(j\) 反而更优。所以随着 \(x\) 增大，最优决策会 向斜率更小的决策移动。

---

三、维护凸壳：弹出队尾的条件

假设队列中已有两条决策 \(a\)（队尾前一个）和 \(b\)（队尾），现在要插入新决策 \(c\)（新来的决策，编号更大，斜率更大）。我们需要判断 \(b\) 是否还会成为最优。

根据上面的分析，\(b\) 和 \(a\) 的优劣分界点为：

\[x_{ab} = -\frac{Y_b - Y_a}{K_b - K_a} \]

\(c\) 和 \(b\) 的分界点为：

\[x_{bc} = -\frac{Y_c - Y_b}{K_c - K_b} \]

3.1 对于 \(\operatorname{opt} = \max\)

• 当 \(x < x_{ab}\) 时，\(a\) 优于 \(b\)；
• 当 \(x \ge x_{ab}\) 时，\(b\) 优于 \(a\)。
• 当 \(x < x_{bc}\) 时，\(b\) 优于 \(c\)；
• 当 \(x \ge x_{bc}\) 时，\(c\) 优于 \(b\)。

如果 \(x_{ab} \ge x_{bc}\)，那么在 \(b\) 刚要优于 \(a\) 的时刻，\(c\) 已经优于 \(b\)，因此 \(b\) 永远不会成为最优，应当弹出。

所以弹出条件是：

\[x_{ab} \ge x_{bc} \]

代入表达式：

\[-\frac{Y_b - Y_a}{K_b - K_a} \ge -\frac{Y_c - Y_b}{K_c - K_b} \]

两边乘以 \(-1\) 并注意分母为正（递增），得：

\[\frac{Y_b - Y_a}{K_b - K_a} \le \frac{Y_c - Y_b}{K_c - K_b} \]

这就是常用的斜率不等式。

交叉相乘（分母均为正）：

\[(Y_b - Y_a)(K_c - K_b) \le (Y_c - Y_b)(K_b - K_a) \]

---

3.2 对于 \(\operatorname{opt} = \min\)

• 当 \(x \le x_{ab}\) 时，\(b\) 优于 \(a\)（因为较小值在左侧）；
• 当 \(x > x_{ab}\) 时，\(a\) 优于 \(b\)。
• 当 \(x \le x_{bc}\) 时，\(c\) 优于 \(b\)；
• 当 \(x > x_{bc}\) 时，\(b\) 优于 \(c\)。

若 \(b\) 要成为最优，需要存在某个区间使 \(b\) 优于 \(a\) 且 \(c\) 优于 \(b\) 不成立。实际上，要求 \(x_{ab} < x_{bc}\)，否则 \(b\) 被夹死。所以弹出条件为：

\[x_{ab} \ge x_{bc} \]

同样得到：

\[\frac{Y_b - Y_a}{K_b - K_a} \ge \frac{Y_c - Y_b}{K_c - K_b} \]

交叉相乘：

\[(Y_b - Y_a)(K_c - K_b) \ge (Y_c - Y_b)(K_b - K_a) \]

---

3.3 统一公式（\(K\) 递增）

操作	\(\operatorname{opt} = \max\)	\(\operatorname{opt} = \min\)
弹出队尾条件	\((Y_b-Y_a)(K_c-K_b) \le (Y_c-Y_b)(K_b-K_a)\)	\((Y_b-Y_a)(K_c-K_b) \ge (Y_c-Y_b)(K_b-K_a)\)
查询弹出队首	比较当前 x 下 q[L] 和 q[L+1] 的值，若后者 \(\ge\) 前者则 L++	若后者 \(\le\) 前者则 L++

四、查询时弹出队首

由于查询点 \(x = X_i\) 单调递增，队首可能不再最优，需要弹出。

• \(\operatorname{opt} = \max\)：比较队首 \(j_1\) 和次队首 \(j_2\) 在 \(x\) 下的函数值，若 \(val(j_2) \ge val(j_1)\)，则 \(j_1\) 被淘汰，指针前移。
• \(\operatorname{opt} = \min\)：若 \(val(j_2) \le val(j_1)\)，则淘汰 \(j_1\)。

---

五、完整推导示例（P3628 特别行动队）

转移方程：

\[dp[i] = \max_{j<i} \left\{ dp[j] + a(S_i-S_j)^2 + b(S_i-S_j) + c \right\} \]

展开：

\[dp[i] = \max_{j<i} \Big\{ \underbrace{dp[j] + aS_j^2 - bS_j}_{Y_j} + \underbrace{(-2aS_j)}_{K_j} \cdot S_i \Big\} + aS_i^2 + bS_i + c \]

这里 \(a<0\)，所以 \(K_j = -2aS_j\) 随 \(j\) 递增（因为 \(-2a>0\)，\(S_j\) 递增）。\(\operatorname{opt} = \max\)。

套用上面的公式，弹出队尾条件为：

\[(Y_b - Y_a)(K_c - K_b) \le (Y_c - Y_b)(K_b - K_a) \]

其中 \(Y_j = dp[j] + aS_j^2 - bS_j\)，\(K_j = -2aS_j\)。

查询时，若 \(val(q[L+1]) \ge val(q[L])\) 则弹出。

---

六、代码模板（以最大值为例）

// 定义 Y[j], K[j]
long long K(int j) { return -2 * a * sum[j]; }
long long Y(int j) { return f[j] + a * sum[j] * sum[j] - b * sum[j]; }

// 计算直线在 x 处的值
long long value(int j, long long x) { return Y(j) + K(j) * x; }

// 维护队列
int head = 1, tail = 0;
q[++tail] = 0;
for (int i = 1; i <= n; ++i) {
    // 弹出队首
    while (head < tail && value(q[head+1], sum[i]) >= value(q[head], sum[i])) head++;
    int j = q[head];
    f[i] = value(j, sum[i]) + a * sum[i] * sum[i] + b * sum[i] + c;

    // 插入新决策，弹出队尾
    while (head < tail) {
        int a = q[tail-1], b = q[tail], c = i;
        __int128 left = (Y(b) - Y(a)) * (K(c) - K(b));
        __int128 right = (Y(c) - Y(b)) * (K(b) - K(a));
        if (left <= right) tail--;
        else break;
    }
    q[++tail] = i;
}

七、通用步骤总结

1. 将转移式化为直线形式
   \(dp_i = \operatorname{opt}\limits_{j<i} (Y_j + K_j X_i) + D_i\)，明确 \(Y_j, K_j, X_i\)。

2. 确定单调性

   • 检查 \(X_i\) 是否单调（通常递增）。
   • 检查 \(K_j\) 是否单调（若递减，可对 \(X\) 取负变为递增，或调整推导中的符号）。

3. 用“k 优于 j”推导弹出条件

   • 写出 \(k\) 优于 \(j\) 的不等式。
   • 解出 \(x\) 的范围，得到分界点表达式。
   • 对三个连续决策 \(a,b,c\)，比较分界点，得出弹出条件（交叉相乘形式）。

4. 查询时弹出队首
   直接比较队首和次队首在当前 \(x\) 下的函数值，按 \(\operatorname{opt}\) 的方向决定是否弹出。

5. 注意溢出
   所有乘法使用 __int128。

---

八、结语

这个方法完全基于代数，无需记忆“上凸壳/下凸壳”的几何方向。每次遇到新题，只需按步骤推导一遍，即可确保方向正确。熟练后，推导过程不超过两分钟，且几乎不会出错。希望这篇笔记能帮助你彻底掌握斜率优化的“方向”难题。

---

P3628 [APIO2010] 特别行动队

#include<bits/stdc++.h>
using namespace std;
const int N=1e6+5;
int n,q[N],L,R,x[N];
#define ll long long
ll f[N],a,b,c,sum[N];

int main(){
	scanf("%d%lld%lld%lld",&n,&a,&b,&c);
	for(int i=1;i<=n;i++){
		scanf("%d",&x[i]);sum[i]=sum[i-1]+x[i];
	}
	L=1;R=1;q[L]=0;
	for(int i=1;i<=n;i++){
		while(L<R&&(-(__int128)2*a*sum[q[L+1]]*sum[i]+f[q[L+1]]+(__int128)a*sum[q[L+1]]*sum[q[L+1]]-b*sum[q[L+1]])>=(-(__int128)2*a*sum[q[L]]*sum[i]+(__int128)a*sum[q[L]]*sum[q[L]]+f[q[L]]-b*sum[q[L]]))L++;
		f[i]=a*sum[i]*sum[i]+b*sum[i]+c+(-2*a*sum[q[L]]*sum[i])+f[q[L]]+a*sum[q[L]]*sum[q[L]]-b*sum[q[L]];
		while(L<R&&__int128((f[q[R]]+(__int128)a*sum[q[R]]*sum[q[R]]-b*sum[q[R]]-f[q[R-1]]-(__int128)a*sum[q[R-1]]*sum[q[R-1]]+b*sum[q[R-1]]))*(-2*a*sum[i]+2*a*sum[q[R]])<=__int128((f[i]+(__int128)a*sum[i]*sum[i]-b*sum[i]-f[q[R]]-(__int128)a*sum[q[R]]*sum[q[R]]+b*sum[q[R]]))*(-2*a*sum[q[R]]+2*a*sum[q[R-1]]))R--;
		q[++R]=i;
	}
	printf("%lld\n",f[n]);
	
	return 0;
}

任务安排1~3

判断时的乘积有可能会超过long long所以需要额外转成__int128比较，但是也不用担心性能问题，因为__int128也没有很慢

#include<stdio.h>
using namespace std;
const int N=3e5+5;
#define ll long long
int n,s,q[N],l,r,t[N],c[N];
ll f[N],ts[N],cs[N];

signed main(){
	scanf("%d%d",&n,&s);
	for(int i=1;i<=n;i++){
		scanf("%d%d",&t[i],&c[i]);
		ts[i]=ts[i-1]+t[i];
		cs[i]=cs[i-1]+c[i];
	}
	q[1]=0;l=1;r=1;
	f[0]=0;
	for(int i=1;i<=n;i++){
		int L=l,R=r,mid;
		while(L<R){
			mid=(L+R)>>1;
			if((f[q[mid]]-f[q[mid+1]])>=(__int128)(s+ts[i])*(cs[q[mid]]-cs[q[mid+1]]))L=mid+1;
			else R=mid;
		}
		f[i]=f[q[L]]-(s+ts[i])*(cs[q[L]])+ts[i]*cs[i]+cs[n]*s;
		while(l<r&&((__int128)(f[q[r]]-f[q[r-1]])*(cs[i]-cs[q[r]])>=(__int128)(f[i]-f[q[r]])*(cs[q[r]]-cs[q[r-1]])))r--;
		q[++r]=i;
	}
	printf("%lld\n",f[n]);
	return 0;
}

P4072 [SDOI2016] 征途

#include<bits/stdc++.h>
using namespace std;
const int N=3e3+5;
int n,m,head,tail,q[N],v[N];
using ll = long long;
const ll INF=1e18;
ll sum[N],f[N][N];

int main(){
	scanf("%d%d",&n,&m);
	for(int i=1;i<=n;i++){
		scanf("%d",&v[i]);
		sum[i]=sum[i-1]+v[i];
	}
	for(int i=0;i<N;i++){f[0][i]=INF;}
	f[0][0]=0;
	for(int d=1;d<=m;d++){
		head=1;q[head]=0;tail=0;
		for(int i=d;i<=n;i++){
			int j=i-1;
			if(f[d-1][j]<INF){
				while(head<tail&&(f[d-1][q[tail]]+(__int128)sum[q[tail]]*sum[q[tail]]-f[d-1][q[tail-1]]-(__int128)sum[q[tail-1]]*sum[q[tail-1]])*(-2*sum[q[tail]]+2*sum[j])>=(f[d-1][j]+(__int128)sum[j]*sum[j]-f[d-1][q[tail]]-(__int128)sum[q[tail]]*sum[q[tail]])*(-2*sum[q[tail-1]]+2*sum[q[tail]]))tail--;
				q[++tail]=j;
			}
			while(head<tail&&f[d-1][q[head]]+(__int128)sum[q[head]]*sum[q[head]]-(__int128)2*sum[q[head]]*sum[i]>=f[d-1][q[head+1]]+(__int128)sum[q[head+1]]*sum[q[head+1]]-(__int128)2*sum[q[head+1]]*sum[i])head++;
			if(head<=tail)f[d][i]=sum[q[head]]*sum[q[head]]-2ll*sum[i]*sum[q[head]]+sum[i]*sum[i]+f[d-1][q[head]];
			else f[d][i]=INF;
		}
	}
	printf("%lld\n",1ll*m*f[m][n]-sum[n]*sum[n]);
	return 0;
}

CF311B Cats Transport

注意不要忘了继承前面的状态

#include<limits.h>
#include<stdio.h>
#include<cstring>
#include<algorithm>
using namespace std;
const int N=1e5+5;
const int P=105;
using ll = long long;
const ll INF=0x3f3f3f3f3f3f3f3f;
int n,m,p,q[N],l,r,d[N],h[N],a[N];
ll f[P][N],sumt[N],sumd[N],t[N];

int main(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=2;i<=n;i++){
		scanf("%d",&d[i]);
		sumd[i]=sumd[i-1]+d[i];
	}
	for(int i=1;i<=m;i++){
		scanf("%d%d",&h[i],&a[i]);
		t[i]=a[i]-sumd[h[i]];
	}
	sort(t+1,t+m+1);
	for(int i=1;i<=m;i++)
		sumt[i]=sumt[i-1]+t[i];
	memset(f,0x3f,sizeof f);
	f[0][0]=0;
	for(int i=1;i<=p;i++){
		for(int j=0;j<=m;j++)
			f[i][j]=f[i-1][j];
		l=0;r=0;q[++l]=0;r++;
		for(int j=1;j<=m;j++){
			while(l<r&&((f[i-1][q[l+1]]+sumt[q[l+1]]-f[i-1][q[l]]-sumt[q[l]])<=1ll*t[j]*(q[l+1]-q[l])))l++;
			f[i][j]=min(f[i][j],f[i-1][q[l]]+sumt[q[l]]-t[j]*q[l]+t[j]*j-sumt[j]);
			if(f[i][j]<INF){
				while(l<r&&((__int128)(f[i-1][q[r]]+sumt[q[r]]-f[i-1][q[r-1]]-sumt[q[r-1]])*(j-q[r])>=(f[i-1][j]+sumt[j]-f[i-1][q[r]]-sumt[q[r]])*(__int128)(q[r]-q[r-1])))r--;
				q[++r]=j;
			}
		}
	}
	printf("%lld\n",f[p][m]);
	
	return 0;
}

Day7 模拟赛

难闷了这次

T1 AT_past201912_n

#include<bits/stdc++.h>
using namespace std;
const int N=1e5+5;
int n,w,c,tot,v[N<<1],L,R,cnt;
struct sg{
	int l,r,p;
}a[N];
using ll = long long;
const ll INF=1e18;
ll cst,dl[N<<1],dr[N<<1],ans=INF;

int main(){
//	freopen("cover.in","r",stdin);
//	freopen("cover.out","w",stdout);
	scanf("%d%d%d",&n,&w,&c);
	for(int i=1;i<=n;i++){
		scanf("%d%d%d",&a[i].l,&a[i].r,&a[i].p);
		v[(i<<1)-1]=a[i].l;v[(i<<1)]=a[i].r;
	}
	n++;v[(n<<1)-1]=0;v[(n<<1)]=w;
	sort(v+1,v+(n<<1)+1);
	tot=unique(v+1,v+(n<<1)+1)-v-1;
	for(int i=1;i<=n;i++){
		a[i].l=lower_bound(v+1,v+tot+1,a[i].l)-v;
		a[i].r=lower_bound(v+1,v+tot+1,a[i].r)-v;
		dr[a[i].l+1]+=a[i].p;dl[a[i].r]-=a[i].p;
	}
	R=0;
	for(L=0;L<=tot;L++){
		cst+=dl[L];
		while(v[R]-v[L]<c&&R<tot){R++;cst+=dr[R];}
		if(v[R]-v[L]>=c){ans=min(ans,cst);}
	}
	printf("%lld\n",ans);
	
	return 0;
}

Day8 高斯消元

基本上是板子的高斯消元

P5027 Barracuda

暴力枚举某个方程冲突，然后求解.我怎么写了这么久？ ？！区区？！

#include<bits/stdc++.h>
using namespace std;
const int N=105;
const double eps=1e-9;
int n,ans;
double a[N][N],b[N],c[N],d[N][N],e[N];

int calc(int n,double a[N][N],double b[N]){
	double ans[N]={0};
	int r=1;
	for(int i=1;i<=n;i++){
		int pos=r;
		for(int j=r+1;j<=n;j++){
			if(fabs(a[pos][i])<fabs(a[j][i]))
				pos=j;
		}
		if(fabs(a[pos][i])<eps){
			continue;
		}
		if(pos!=r){
			swap(a[pos],a[r]);
			swap(b[pos],b[r]);
		}
		double div=a[r][i];
		for(int k=i;k<=n;k++) a[r][k]/=div;
		b[r]/=div;
		for(int j=1;j<=n;j++){
			if(j==r) continue;
			double rate=a[j][i];
			for(int k=i;k<=n;k++) a[j][k]-=a[r][k]*rate;
			b[j]-=b[r]*rate;
		}
		++r;
	}
	bool f=0;
	for(int i=r;i<=n;i++){
		bool zl=1;
		for(int j=1;j<=n;j++){
			if(fabs(a[i][j])>eps){
				zl=0;break;
			}
		}
		if(zl&&fabs(b[i])>eps){
			return -1;
		}
		if(zl)f=1;
	}
	if(f||r-1<n){
		return -1;
	}
	for(int i=1;i<=n;i++) ans[i]=b[i];
	double tmp=0.0;
	int res=0,cnt=0;
	for(int i=1;i<=n;i++){
		if(ans[i]<1-eps)return -1;
		if(fabs(ans[i]-round(ans[i]))>eps)return -1;
		ans[i]=round(ans[i]);
		if(fabs(ans[i]-tmp)<=eps){++cnt;}
		if(ans[i]>tmp){
			tmp=ans[i];res=i;cnt=1;
		}
	}
	if(cnt>1)return -1;
	return res;
}

int main(){
	scanf("%d",&n);
	for(int i=1;i<=n+1;i++){
		int m;scanf("%d",&m);
		for(int j=1;j<=m;j++){
			int x;scanf("%d",&x);
			a[i][x]+=1;
		}
		scanf("%lf",&b[i]);
	}
	for(int i=1;i<=n+1;i++){
		int t=0;memset(d,0,sizeof d);memset(e,0,sizeof e);
		for(int j=1;j<=n+1;j++){
			if(j==i)continue;
			++t;
			for(int k=1;k<=n;k++){
				d[t][k]=a[j][k];
			}
			e[t]=b[j];
		}
		int tmp=calc(n,d,e);
		if(tmp>0&&ans!=0){
			printf("illegal\n");
			return 0;
		}
		else if(tmp>0)ans=tmp;
	}
	if(ans)printf("%d\n",ans);
	else printf("illegal\n");
	
	return 0;
}

P2455 [SDOI2006] 线性方程组

模板题，但是要求输出无解和无穷解的情况，很坑，需要厘清代码实现逻辑，我就没有所以WA了很久很不错的一道题，使人理解深刻.但是这里我他妈的写了很久一直90分最后使用矩阵的秩以及在循环外面单独在循环判断才可以过

矩阵的秩的通俗解释:简单说，一个矩阵的秩，就是它里面真正有效、独立的方程个数.

不用管形式化定义至少目前来说是这样

#include<bits/stdc++.h>
using namespace std;
const int N=105;
const double eps=1e-12;
int n;
double b[N],a[N][N],ans[N];

int main(){
	scanf("%d",&n);
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++){
			scanf("%lf",&a[i][j]);
		}
		scanf("%lf",&b[i]);
	}
	int r=1;//r即矩阵的秩
	for(int i=1;i<=n;i++){
		int pos=r;
		for(int j=r;j<=n;j++){
			if(fabs(a[pos][i])-fabs(a[j][i])<eps)
				pos=j;
		}
		if(fabs(a[pos][i])<eps){
			continue;
		}
		if(pos!=r){
			swap(a[pos],a[r]);
			swap(b[pos],b[r]);
		}
		for(int j=r+1;j<=n;j++){
			double rate=a[j][i]/a[r][i];
			for(int k=i;k<=n;k++)
				a[j][k]-=a[r][k]*rate;
			b[j]-=b[r]*rate;
		}
		++r;
	}
	bool f=0;
	for(int i=r;i<=n;i++){
		bool zl=1;
		for(int j=1;j<=n;j++){
			if(fabs(a[i][j])>eps){
				zl=0;break;
			}
		}
		if(zl&&fabs(b[i])>eps){
			printf("-1\n");return 0;
		}
		if(zl)f=1;
	}
	if(f||r-1<n){
		printf("0\n");return 0;
	}
	ans[n]=b[n]/a[n][n];
	for(int i=n-1;i>=1;i--){
		double s=b[i];
		for(int j=i+1;j<=n;j++)
			s-=a[i][j]*ans[j];
		ans[i]=s/a[i][i];
	}
	for(int i=1;i<=n;i++)printf("x%d=%.3lf\n",i,fabs(ans[i])<eps? 0.0:ans[i]);//防止输出"-0.000"这样的答案(虽然说有spj就不用管这个但我还是写了,原因是测样例时看起来也太奇怪了)
	return 0;
}

高斯消元求动态规划式子的值

总感觉这类题多的一批，每次还必是 \(\color{purple}{省选}\) 难度

P2973 [USACO10HOL] Driving Out the Piggies G

感觉和随机游走挺像的，对于本题，求被摧毁的概率，由于每个炸弹只有一次爆炸机会，所以我们可以轻松得到，炸弹在 \(i\) 处爆炸的概率就等于其经过 \(i\) 的期望次数 \(\times\) 爆炸概率，为什么呢？

这里引用来自洛谷题解区的证明:

其它题解提到，在一个点爆炸的概率 = 这个点的期望经过次数 × 每次炸弹爆炸的概率（即 \(\frac{p}{q}\)）。

这样感性理解一下可能能觉得是对的。这里写一下一个简单的，不那么偏感性的解释：

显然炸弹爆炸的次数的取值只有 \(\mathcal{\text{{0,1}}}\)。这样我们发现，炸弹的期望爆炸次数与爆炸概率 在数值上相等。从而

在一个点爆炸的概率 = 期望在这个点爆炸的次数 = 期望经过这个点且在这个点爆炸的次数

当随机变量 \(\mathcal{X}\) 和 \(\mathcal{Y}\) 互相独立时，\(\mathcal{E(XY)=E(X)\times E(Y)}\) 。定义 \(\mathcal{X}\) 为经过这个点的次数，\(\mathcal{Y}\) 为每次爆炸次数（显然同样取值 \(\mathcal{\text{{0,1}}}\) ），则所求即

\(\mathcal{E(XY)=E(X)\times E(Y)}\) 同样，由 \(\mathcal{Y} \in \mathcal{\text{{0,1}}}\) 知，\(\mathcal{E(Y)}\) 数值上等于每次爆炸概率，即 \(\frac{p}{q}\) 。

所以，在一个点爆炸的概率，数值上等于在这个点的期望爆炸次数 \(\mathcal{E(XY)}\) ，数值上等于这个点的期望经过次数 \(\mathcal{E(X)}\) 乘每次爆炸概率 \(\frac{p}{q}\)

所以我们定义 \(g_i\) 表示从 \(1\) 出发经过点 \(i\) 的期望次数,显然有:

\[\mathcal{g_i} = [i=1]+\sum\limits_{(i,v)\in E} g_v ·\frac{1-P}{du_v} \]

其中 \(P=\frac{p}{q}\)，\(du_v\) 即为节点 \(v\) 的度数, \(g_1\) 要加一是因为我们是从 \(1\) 号节点开始，所以它开始就被访问了一次.

而我们最后要求的答案就等于 \(P·g_i\). 时间复杂度 \(\mathcal{O(n^4)}\) 正好通过本题.

#include<bits/stdc++.h>
using namespace std;
const int N=305;
const double eps=1e-6;
int n,m,du[N];
double p,q,P,a[N][N],b[N],ans[N];
vector<int> g[N];

void gauss(){
	int r=1;
	for(int i=1;i<=n;i++){
		int pos=r;
		for(int j=r+1;j<=n;j++){
			if(fabs(a[pos][i])-fabs(a[j][i])<eps)pos=j;
		}
		if(pos!=r){
			swap(a[pos],a[r]);
			swap(b[pos],b[r]);
		}
		for(int j=r+1;j<=n;j++){
			double rate=a[j][i]/a[r][i];
			for(int k=i;k<=n;k++)
				a[j][k]-=a[r][k]*rate;
			b[j]-=b[r]*rate;
		}
		++r;
	}
	ans[n]=b[n]/a[n][n];
	for(int i=n-1;i>=1;i--){
		double s=b[i];
		for(int j=i+1;j<=n;j++)
			s-=a[i][j]*ans[j];
		ans[i]=s/a[i][i];
	}
}

int main(){
	scanf("%d%d%lf%lf",&n,&m,&p,&q);P=(p/q);
	for(int i=1;i<=m;i++){
		int u,v;scanf("%d%d",&u,&v);du[u]++;du[v]++;
		g[u].push_back(v);g[v].push_back(u);
	}
	for(int u=1;u<=n;u++){
		a[u][u]=1.0;
		for(auto v:g[u]){
			a[u][v]=-(1-P)/du[v];
		}
		if(u==1)b[u]=P;
	}
	gauss();
	for(int i=1;i<=n;i++)
		printf("%.10lf\n",ans[i]);
	
	return 0;
}

[HNOI2013] 游走

#include<bits/stdc++.h>
using namespace std;
const int N=505;
const int M=1.25e5+5;
const double eps=1e-12;
int n,m,st[M],ed[M],head[N],tot,du[N];
struct Edge{
	int nxt,to;
}edge[M<<1];
double a[N][N],f[N],g[M],ans;
void add(int u,int v){
	edge[++tot].to=v;
	edge[tot].nxt=head[u];
	head[u]=tot;
}

inline void gauss(int n){
    for(int i=1;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(fabs(a[r][i])<eps)continue;
        if(r!=i)
			for(int j=1;j<=n+1;j++)
				swap(a[i][j],a[r][j]);
        for(int j=1;j<=n;j++){
            if(j!=i&&fabs(a[j][i])>eps){
                double rate=a[j][i]/a[i][i];
                for(int k=i;k<=n+1;k++)
					a[j][k]-=a[i][k]*rate;
            }
        }
    }
    for(int i=1;i<=n;i++)
		f[i]=(fabs(a[i][i])>eps)?a[i][n+1]/a[i][i]:0;
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=1;i<=m;i++){
        scanf("%d%d",&st[i],&ed[i]);
        add(st[i],ed[i]);add(ed[i],st[i]);
        du[st[i]]++;du[ed[i]]++;
    }
    for(int i=1;i<n;i++){
        a[i][i]=1.0;
        for(int j=head[i];j;j=edge[j].nxt){
            int v=edge[j].to;
            if(v!=n)a[i][v]-=1.0/du[v];
        }
        a[i][n]=(i==1)?1.0:0.0;
    }
    gauss(n-1);
    for(int i=1;i<=m;i++){
        int u=st[i],v=ed[i];
        if(u!=n)g[i]+=f[u]/du[u];
        if(v!=n)g[i]+=f[v]/du[v];
    }
    sort(g+1,g+m+1);
    for(int i=1;i<=m;i++)ans+=g[i]*(m-i+1);
    printf("%.3lf\n",ans);
    return 0;
}

矩阵快速幂加速递推

P3913 [HNOI2013] GT考试

先考虑普通动态规划。设状态j表示当前已经构造出的数字串的最长后缀，恰好等于不吉利串的长度为 \(j\) 的前缀，其中 \(0≤j<M\) 。状态不能到达 \(M\)，因为那意味着已经出现了完整的不吉利串。
对每个状态 \(j\) 和每个新加入的数字 \(c\)，需要求加入后新的最长匹配前缀长度。这个过程和KMP自动机一致：如果当前字符不能继续匹配，就沿着前缀函数回退，直到可以匹配或回到 \(0\)。若加入后长度变成 \(M\)，这条转移被丢弃；否则把对应转移计数加一。
这样可以得到一个 \(M×M\) 的转移矩阵 \(A\)，其中 \(A_{i,j}\) 表示从状态 \(i\) 加入一个数字后转到状态 \(j\) 的数字选择数。并且，和前面一样 \(dp_{i,j}\) 表示前 \(i\) 个字符的后缀有 \(j\) 与不吉利串的前缀匹配的方案数，而且我们仍然保证 \(0\le j < M\).初始时处于状态 \(0\) 最多只能和不吉利串的前 \(0\) 个字符匹配一次,即两个空串匹配有 \(1\) 的方案数，所以记行向量 \(dp_0 = [1,0,...,0]\).
那么构造 \(i\) 位数字串后，根据矩阵乘法，我们有转移方程：

\[dp_i = dp_{i-1}\times A. \]

即

\[dp_i = dp_0\times A^i. \]

答案为即 \(dp_N\) 中所有状态的方案数之和。
由于 \(N\) 很大，使用矩阵快速幂求 \(A_N\)。本题 \(M≤20\)，矩阵规模很小，可以直接三重循环乘法.时间复杂度 \(\mathcal{O}(M^3\log N)\). 可通过本题,实现时注意数组下标的含义要对应.应该没人和我一样下标错了好几次吧

#include<bits/stdc++.h>
using namespace std;
const int M=25;
int n,m,mod,f0[M],p[M],fn[M];
char ch[M];
struct matrix{
	int a[M][M]={0};
	matrix operator*(const matrix& b) const{
		matrix c;
		for(int i=0;i<m;i++){
			for(int k=0;k<m;k++){
				if(a[i][k]==0)continue;
				for(int j=0;j<m;j++){
					c.a[i][j]=(c.a[i][j]+a[i][k]*b.a[k][j]%mod)%mod;
				}
			}
		}
		return c;
	}
}A;

matrix fastpow(matrix a,int b){
	matrix res;
	if(b){res=a;b--;}
	for(;b;b>>=1){
		if(b&1)res=res*a;
		a=a*a;
	}
	return res;
}

int main(){
	scanf("%d%d%d",&n,&m,&mod);
	scanf("%s",ch);
	for(int i=1,j=0;i<m;i++){
		while(j&&ch[i]!=ch[j])j=p[j-1];
		if(ch[j]==ch[i])++j;
		p[i]=j;
	}
	for(int i=0;i<m;i++){
		for(int c='0';c<='9';c++){
			int j=i;
			while(j&&ch[j]!=c)j=p[j-1];
			if(ch[j]==c)++j;
			if(j<m)A.a[i][j]=(A.a[i][j]+1)%mod;
		}
	}
	A=fastpow(A,n);
	f0[0]=1;
	for(int j=0;j<m;j++){
		for(int i=0;i<m;i++){
			fn[j]=(fn[j]+f0[i]*A.a[i][j])%mod;
		}
	}
	int ans=0;
	for(int i=0;i<m;i++)
		ans=(ans+fn[i])%mod;
	printf("%d\n",ans);
	
	return 0;
}