隐式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）