隐式dijkstra:在状态集合中用优先队列求前k小(待续)
这种技巧是挺久以前接触的了,最近又突然遇到几道新题,于是总结了一下体会。
要注意,这种算法适用的前提是,标题所述的“状态集合”的size大到不可枚举(否则直接枚举就行了qaq),而且k一般是在1e6这个数量级以下。
前置技能:Dijkstra算法,及其思想和正确性证明。
传送门1:思想和正确性证明
传送门2:优先队列优化dijkstra
先看一个问题:
给1<m<10个长度为n<1e5的整数序列,从每个序列选一个数相加,求所得的和中第k<1e5大的。
(首先显而易见要把每个序列排序,从大到小)
考虑m=2的情况:
二分答案A,就可以对序列1中的每一个元素,找到能够使总和在A以上的、序列2中可以与它配对的元素集合。易知这个集合是序列2的一个前缀,且前缀的长度随i递增,故对一个A求出所有这样的前缀只需要线性时间。而总和大于A的方案数,就是这些前缀的长度之和。
上代码:
#include <bits/stdc++.h> using namespace std; #define iinf 1000000000 #define linf 1000000000000000000LL #define ulinf 10000000000000000000ull #define MOD1 1000000007LL #define mpr make_pair typedef long long LL; typedef unsigned long long ULL; typedef unsigned long UL; typedef unsigned short US; typedef pair < int , int > pii; clock_t __stt; inline void TStart(){__stt=clock();} inline void TReport(){printf("\nTaken Time : %.3lf sec\n",(double)(clock()-__stt)/CLOCKS_PER_SEC);} template < typename T > T MIN(T a,T b){return a<b?a:b;} template < typename T > T MAX(T a,T b){return a>b?a:b;} template < typename T > T ABS(T a){return a>0?a:(-a);} template < typename T > void UMIN(T &a,T b){if(b<a) a=b;} template < typename T > void UMAX(T &a,T b){if(b>a) a=b;} int n,m,k,a[2][100005]; bool check(int v){ int i,j,t,cnt; for(t=0;t<n && a[0][t]+a[1][0]<v;++t); cnt=n-t; for(i=1;i<n;++i){ for(;t<n && a[0][t]+a[1][i]<v;++t); cnt+=n-t; if(cnt>=k) return 1; } return cnt>=k; } int main(){ // inputting start // 数据结构记得初始化! n,m别写反! scanf("%d%d%d",&n,&m,&k); //m=2 int i,j; for(i=0;i<2;++i){ for(j=0;j<n;++j){ scanf("%d",&a[i][j]); } } #ifdef LOCAL TStart(); #endif // calculation start // 数据结构记得初始化! n,m别写反! sort(a[0],a[0]+n); sort(a[1],a[1]+n); reverse(a[1],a[1]+n); int low=0,high=iinf,mid; while(low<high){ mid=((low+high+1)>>1); if(check(mid)) low=mid; else high=mid-1; } printf("%d\n",low); #ifdef LOCAL TReport(); #endif return 0; }
这种方法不是本文的重点,但是很大一部分可以用马上会介绍的隐式dijkstra解决的问题,都可以用这种二分的方法,以一个略差的性能解决。
有的问题看似可以用隐式dijkstra,其实却会出大问题。碰到这种情况,可以考虑改用这类二分答案的做法。
推荐下面这道题。
USACO 2016DEC PLATINUM T3
题解
对于一般情况(m<=10),就要换方法了。
(敲黑板!重点重点!)
我们在脑海中建一张图,每个结点对应着一种选取方案。方案A到B有一条有向边,当且仅当A对应的总和大于B的总和,边权是两者总和的差的绝对值。
记所有序列最大数组成S状态,那么对任意T状态:cost(T)=cost(S)-dist(S,T)
那么我们要求的,就是到S状态距离第k短的点(包括S自己)。注意到虽然这是个DAG,但因为结点过多无法DP。
所以,考虑dijkstra。
先上结论:
dijkstra在正权图上运行时,优先队列每次弹出的结点到S的距离,一定是递增的。
证明:
若有dist(S,P)>dist(S,Q)而P先于Q弹出,则:
在P弹出前的瞬间,因为P是优先队列中距S最近的,所以Q到S比优先队列中任意状态到S更近,且Q在P弹出前不在优先队列里,也从未被压入过。
所以Q不可能由优先队列中的状态经过若干松弛操作而得到。
故Q不可能在P弹出后被压入队列,也就不可能在P之后弹出。
产生矛盾,证毕。
推论:
对任意的k,优先队列里最先弹出的k个结点,一定是到源点最近的k个点。
我会做啦!啊哈哈哈哈哈!
跑一遍dijkstra,优先队列弹出的第k个点就是答案!
且慢,算一下复杂度。
记各结点的平均度数为d,因为做了k轮松弛,每轮压入了d个新结点,故总时间复杂度为:
O(kd*log(k)),约等于1e56。emmmm……
(敲黑板!重点又来了!)
这种算法,优化的思路主要有两种,在此先用第一种:
优化连边!哇哈哈哈哈哈!
易知,对于固定的T,任意S-T路径的权值和全都相等。
所以,如果删掉一些边,使得图的连通性不变,那么答案也不变。
换句话说,要删掉一些边,使得S到每个点仍有至少一条路径。
所以修改连边策略:
A到B有一条有向边,仅当A对应的总和大于B的总和,边权是两者总和的差的绝对值。
此外,必须满足,A和B对应的方案,只在一个序列里选的数下标不同(其它m-1个序列里选的都完全相同),而且,所选的两个下标不同的数,也一定是相邻元素。(数组已排序)
要删掉一些边,使得S到每个点仍有至少一条路径。
满足了没?满足了。
而现在,d<=m。所以总复杂度降为O(mk*log(k)),bingo!
具体实现的hint:
实际上,完全可以枚举得到每一个点的所有邻居,于是就没有必要建立邻接表了。
(我不知道这种算法的正式名称(可能并不存在?),所以就在本文中叫它“隐式dijkstra”了。)
上代码:
#include <bits/stdc++.h> using namespace std; #define iinf 1000000000 #define linf 1000000000000000000LL #define ulinf 10000000000000000000ull #define MOD1 1000000007LL #define mpr make_pair typedef long long LL; typedef unsigned long long ULL; typedef unsigned long UL; typedef unsigned short US; typedef pair < int , int > pii; clock_t __stt; inline void TStart(){__stt=clock();} inline void TReport(){printf("\nTaken Time : %.3lf sec\n",(double)(clock()-__stt)/CLOCKS_PER_SEC);} template < typename T > T MIN(T a,T b){return a<b?a:b;} template < typename T > T MAX(T a,T b){return a>b?a:b;} template < typename T > T ABS(T a){return a>0?a:(-a);} template < typename T > void UMIN(T &a,T b){if(b<a) a=b;} template < typename T > void UMAX(T &a,T b){if(b>a) a=b;} int n,m,k,a[2][100005]; struct P{ int x,y; int val(){return a[0][x]+a[1][y];} bool operator <(P b) const{ return a[0][x]+a[1][y]<b.val(); } }; P make_P(int x,int y){ P R; R.x=x;R.y=y; return R; } priority_queue < P > pq; map < pii , int > vis; int main(){ // inputting start // 数据结构记得初始化! n,m别写反! scanf("%d%d%d",&n,&m,&k); //m=2 int i,j; for(i=0;i<2;++i){ for(j=0;j<n;++j){ scanf("%d",&a[i][j]); } } #ifdef LOCAL TStart(); #endif // calculation start // 数据结构记得初始化! n,m别写反! sort(a[0],a[0]+n); reverse(a[0],a[0]+n); sort(a[1],a[1]+n); reverse(a[1],a[1]+n); pq.push(make_P(0,0)); while(!pq.empty()){ P cur=pq.top(); pq.pop(); --k; if(!k){ printf("%d\n",cur.val()); return 0; } if(cur.x<n-1 && !vis[mpr(cur.x+1,cur.y)]){ vis[mpr(cur.x+1,cur.y)]=1; pq.push(make_P(cur.x+1,cur.y)); } if(cur.y<n-1 && !vis[mpr(cur.x,cur.y+1)]){ vis[mpr(cur.x,cur.y+1)]=1; pq.push(make_P(cur.x,cur.y+1)); } } #ifdef LOCAL TReport(); #endif return 0; }
小结1:隐式dijkstra的适用条件
1、合理的数据范围
要注意,这种算法适用的前提是,标题所述的“状态集合”的size大到不可枚举(否则直接枚举就行了qaq),而且k一般是在1e6这个数量级以下。
这里的“状态集合”指的是构出的图中,所有结点的集合。
2、易于表示、比较的状态
优先队列里的操作是基于比较的。
如果比较大小的复杂度过高会TLE,如果存储状态的空间复杂度过大会MLE。
3、正权图
边权必须都是非负数。
实践出真知:SGU421
题意:
给长度为n(<1e4)的整数(正负均可)数列,选m(<13)个数相乘,问第k大的乘积。
题解:
先把数列按正负分成两个数组,再分别排序。
以下是隐式dijkstra的几个要素:
状态含义:
一个状态代表一种选取方案,即一个正数子集和一个负数子集的二元组。
状态间的大小关系:
即乘积的大小关系。先比较符号(负数数量的奇偶性),再比绝对值(高精度)。
起始状态 S:
考虑使用多个起始状态。每个S是负数集合大小为x的状态中最大的。x为任意自然数。
连边方案:
如果状态u,v只有一个选的数不同,而这两个不同的数是相邻元素,则从较大状态向较小状态连边,边权是权值相除的商。
答案就是优先队列弹出的第k个状态的值。
正确性证明:
注意到这张图不是弱连通的,但是每个弱连通分量(包含的状态拥有相同的负数子集大小)都恰有一个起始状态,故正确性依然保持。
上代码:
#include <bits/stdc++.h> using namespace std; #define iinf 2000000000 #define linf 1000000000000000000LL #define ulinf 10000000000000000000ull #define MOD1 1000000007LL #define mpr make_pair typedef long long LL; typedef unsigned long long ULL; typedef unsigned long UL; typedef unsigned int US; typedef pair < int , int > pii; clock_t __stt; inline void TStart(){__stt=clock();} inline void TReport(){printf("\nTaken Time : %.3lf sec\n",(double)(clock()-__stt)/CLOCKS_PER_SEC);} template < typename T > T MIN(T a,T b){return a<b?a:b;} template < typename T > T MAX(T a,T b){return a>b?a:b;} template < typename T > T ABS(T a){return a>0?a:(-a);} template < typename T > void UMIN(T &a,T b){if(b<a) a=b;} template < typename T > void UMAX(T &a,T b){if(b>a) a=b;} int n,m,k; vector < int > neg,pos; struct bigint{ int len,cnt0; int d[85]; void init(){ memset(d,0,sizeof(d)); len=1; d[0]=1; cnt0=0; } void reduct(){ while(len>1 && !d[len-1]) --len; } void multiply(int v){ if(!v){ ++cnt0; return; } int i; for(i=0;i<len;++i) d[i]*=v; for(i=0;i<len;++i){ d[i+1]+=d[i]/10; d[i]%=10; } while(d[len]){ d[len+1]+=d[len]/10; d[len++]%=10; } } void divide(int v){ if(!v){ --cnt0; return; } int i,j,c=0; for(i=len-1;i>=0;--i){ c=c*10+d[i]; d[i]=0; if(c>=v){ d[i]=c/v; c%=v; } } for(i=0;i<len;++i){ d[i+1]+=d[i]/10; d[i]%=10; } while(d[len]){ d[len+1]+=d[len]/10; d[len++]%=10; } reduct(); } void print(bool sig){ if(cnt0){ printf("0\n"); return; } int i; reduct(); if(sig && !(len==1&&d[0]==0)) printf("-"); for(i=len-1;i>=0;--i) printf("%d",d[i]); printf("\n"); } }; bool operator <(bigint &A,bigint &B){ if(A.cnt0) return !B.cnt0; if(B.cnt0) return 0; if(A.len!=B.len) return A.len<B.len; int i; for(i=A.len-1;i>=0;--i){ if(A.d[i]!=B.d[i]) return A.d[i]<B.d[i]; } return 0; } struct state{ int vp[15],vn[15],cp,cn; bigint val; state(){ cp=cn=0; val.init(); } bool sign(){ return cn&1; } bool editp(int p,int d){ if(vp[p]+d<0 || vp[p]+d>=(int)pos.size()) return 0; if((p && vp[p-1]==vp[p]+d)||(p<cp-1 && vp[p+1]==vp[p]+d)) return 0; val.divide(pos[vp[p]]); vp[p]+=d; val.multiply(pos[vp[p]]); return 1; } bool editn(int p,int d){ if(vn[p]+d<0 || vn[p]+d>=(int)neg.size()) return 0; if((p && vn[p-1]==vn[p]+d)||(p<cn-1 && vn[p+1]==vn[p]+d)) return 0; val.divide(neg[vn[p]]); vn[p]+=d; val.multiply(neg[vn[p]]); return 1; } int super_cdd(){ int ret=cp*101+cn,i; for(i=0;i<cp;++i){ ret=ret*101+vp[i]; } for(i=0;i<cn;++i){ ret=ret*101+vn[i]; } return ret; } }; const bool operator <(state A,state B){ if(A.sign()!=B.sign()) return A.sign()>B.sign(); A.val.reduct();B.val.reduct(); return (A.sign()?(B.val<A.val):(A.val<B.val)); } priority_queue < state > pq; map < int , bool > vis; int main(){ // inputting start // 数据结构记得初始化! n,m别写反! scanf("%d%d%d",&n,&m,&k); int i,j; for(i=0;i<n;++i){ scanf("%d",&j); if(j<0){ neg.push_back(-j); } else{ pos.push_back(j); } } #ifdef LOCAL TStart(); #endif // calculation start // 数据结构记得初始化! n,m别写反! sort(pos.begin(),pos.end()); sort(neg.begin(),neg.end()); for(i=0;i<=m;i+=2){ state tmp; for(j=0;j<i && j<(int)neg.size();++j){ tmp.vn[tmp.cn++]=(int)neg.size()-1-j; tmp.val.multiply(neg[(int)neg.size()-1-j]); } if(j>=i){ for(j=0;j<m-i && j<(int)pos.size();++j){ tmp.vp[tmp.cp++]=(int)pos.size()-1-j; tmp.val.multiply(pos[(int)pos.size()-1-j]); } if(j>=m-i){ vis[tmp.super_cdd()]=1; pq.push(tmp); } } } for(i=1;i<=m;i+=2){ state tmp; for(j=0;j<i && j<(int)neg.size();++j){ tmp.vn[tmp.cn++]=j; tmp.val.multiply(neg[j]); } if(j>=i){ for(j=0;j<m-i && j<(int)pos.size();++j){ tmp.vp[tmp.cp++]=j; tmp.val.multiply(pos[j]); } if(j>=m-i){ vis[tmp.super_cdd()]=1; pq.push(tmp); } } } while(k--){ state cur=pq.top(); pq.pop(); if(!k){ cur.val.print(cur.sign()); break; } for(i=0;i<cur.cp;++i){ if(cur.editp(i,(cur.sign()?1:-1))){ if(!vis[cur.super_cdd()]){ vis[cur.super_cdd()]=1; pq.push(cur); } cur.editp(i,(cur.sign()?-1:1)); } } for(i=0;i<cur.cn;++i){ if(cur.editn(i,(cur.sign()?1:-1))){ if(!vis[cur.super_cdd()]){ vis[cur.super_cdd()]=1; pq.push(cur); } cur.editn(i,(cur.sign()?-1:1)); } } } #ifdef LOCAL TReport(); #endif return 0; }
实践出真知2:待更
这题是我打算出出来的一个idea……先等我把它放到oj上再说吧qaq……
(其实我觉得这可能是本文里最妙的一道题了qaq)