统计的力量——线段树全接触

by zkw, Jul 6, 2011, 5:40 PM

　　主要内容是: 一种快速, 易编写的线段树, 标记永久化的思想, 线段树在统计上的应用, 使用条件, 以及部分取代其他数据结构如平衡树的讨论.

　　下载幻灯片

This post has been edited 3 times. Last edited by zkw, Jun 21, 2014, 4:21 AM

BST 拓展与伸展树 (Splay) 一日通

by zkw, Jul 6, 2011, 5:34 PM

　　BST (二分检索树) 是 OI 中一个常用的数据结构, 主要支持的操作是动态维护一个有序表, 从而支持字典, 前驱, 后继, 中序遍历, 优先队列等等多种操作. 如果在域中维护了子树的 size, 还可以支持查找第 k 大的数, 询问名次等等附加功能. 伸展树 (Sleator and Tarjan, 1985) 是 BST 的一个变种, 可以自调整以维护平衡, 并支持根据需要随时把任意节点旋转到根 (称为 splay 操作), 从而很好的支持了分裂合并等操作, 从某种意义上 (详见下文) 也简化了思维和编程的复杂度. 由于 splay 操作的迅速而优美, 伸展树不仅可以维护有序表, 还可以放弃关键字的有序性, 像块状链表一样维护一个一般的序列 (这是一般的 BST 甚至其他平衡二叉树难于做到的), 支持块状链表的大多数方法, 并拥有更低的理论复杂度和实际代码量. 另外, splay 是唯一支持稳定排序的平衡二叉树, 可以妥善处理关键字相同的元素 (按照插入顺序有序), 在输出时可以选择返回其中最近插入、最远插入的, 或是直接返回整个区间.

　　1. 伸展树的基本操作与应用 点击阅读

　　2. 自顶向下伸展树 点击阅读

　　3. Code Trick 点击阅读

　　4. 复杂度 点击阅读

　　5. 应用 点击阅读

　　6. C++ 实现 点击查看程序代码

#include <cstdio>
using namespace std;
const int maxint=~0U>>1;

struct node
{
    int key,size;
    node *c[2];
    node():key(0),size(0){c[0]=c[1]=this;}
    node(int key_,node* c0_,node* c1_):
        key(key_){c[0]=c0_;c[1]=c1_;}
    node* rz(){return size=c[0]->size+c[1]->size+1,this;}
} Tnull,*null=&Tnull;

struct splay
{
    node *root;
    splay()
    {
        root=(new node(*null))->rz();
        root->key=maxint;
    }
    void zig(bool d)
    {
        node *t=root->c[d];
        root->c[d]=null->c[d];
        null->c[d]=root;
        root=t;
    }
    void zigzig(bool d)
    {
        node *t=root->c[d]->c[d];
        root->c[d]->c[d]=null->c[d];
        null->c[d]=root->c[d];
        root->c[d]=null->c[d]->c[!d];
        null->c[d]->c[!d]=root->rz();
        root=t;
    }
    void finish(bool d)
    {
        node *t=null->c[d],*p=root->c[!d];
        while(t!=null)
        {
            t=null->c[d]->c[d];
            null->c[d]->c[d]=p;
            p=null->c[d]->rz();
            null->c[d]=t;
        }
        root->c[!d]=p;
    }
    void select(int k)
    {
        int t;
        for(;;)
        {
            bool d=k>(t=root->c[0]->size);
            if(k==t||root->c[d]==null)break;
            if(d)k-=t+1;
            bool dd=k>(t=root->c[d]->c[0]->size);
            if(k==t||root->c[d]->c[dd]==null){zig(d);break;}
            if(dd)k-=t+1;
            d!=dd?zig(d),zig(dd):zigzig(d);
        }
        finish(0),finish(1);
        root->rz();
    }
    void search(int x)
    {
        for(;;)
        {
            bool d=x>root->key;
            if(root->c[d]==null)break;
            bool dd=x>root->c[d]->key;
            if(root->c[d]->c[dd]==null){zig(d);break;}
            d!=dd?zig(d),zig(dd):zigzig(d);
        }
        finish(0),finish(1);
        root->rz();
        if(x>root->key)select(root->c[0]->size+1);
    }
    void ins(int x)
    {
        search(x);
        node *oldroot=root;
        root=new node(x,oldroot->c[0],oldroot);
        oldroot->c[0]=null;
        oldroot->rz();
        root->rz();
    }
    void del(int x)
    {
        search(x);
        node *oldroot=root;
        root=root->c[1];
        select(0);
        root->c[0]=oldroot->c[0];
        root->rz();
        delete oldroot;
    }
    int sel(int k){return select(k-1),root->key;}
    int ran(int x){return search(x),root->c[0]->size+1;}
} sp;

int main()
{
    for(;;)
    {
        char cmd;
        int num;
        scanf(" %c%d",&cmd,&num);
        switch(cmd)
        {
            case'i':sp.ins(num);break;
            case'd':sp.del(num);break;
            case's':printf("%d\n",sp.sel(num));break;
            case'r':printf("%d\n",sp.ran(num));break;
        }
    }
}

This post has been edited 1 time. Last edited by zkw, Jul 6, 2011, 5:36 PM

SAP (Shortest Augmenting Paths) 最大流算法例程

by zkw, Jul 6, 2011, 5:17 PM

　　偶然翻到还在用 Pascal 时写的 SAP 心得, 发现还有很多人关注, 既然现在我已经在用 C++ 了, 那么这里给出一个C++的版本, 已经加上了当前弧, GAP, 多路增广. 毕竟是C++, 程序比Pascal写出来要短一些.

　　点击查看程序代码

#include <cstdio>
#include <cstring>
using namespace std;
const int V=220,E=220;

int n,m,h[V],vh[V];

struct etype
{
    int t,u;
    etype *next,*pair;
    etype(){}
    etype(int t_,int u_,etype* next_):t(t_),u(u_),next(next_){}
    void* operator new(unsigned,void* p){return p;}
} *e[V],*eb[V],Te[E+E],*Pe=Te;

int aug(int no,int m)
{
    if(no==n)return m;
    int l=m;
    for(etype *&i=e[no];i;i=i->next)
        if(i->u && h[i->t]+1==h[no])
        {
            int d=aug(i->t,l<i->u?l:i->u);
            i->u-=d,i->pair->u+=d,l-=d;
            if(h[1]==n || !l)return m-l;
        }
    int minh=n;
    for(etype *i=e[no]=eb[no];i;i=i->next)if(i->u)
        if(h[i->t]+1<minh)minh=h[i->t]+1;
    if(!--vh[h[no]])h[1]=n;else ++vh[h[no]=minh];
    return m-l;
}
    
int main()
{
    freopen("ditch.in","r",stdin);
    freopen("ditch.out","w",stdout);
    scanf("%d %d",&m,&n);
    memset(e,0,sizeof(e));
    while(m--)
    {
        int s,t,u;
        scanf("%d %d %d",&s,&t,&u);
        e[s]=new(Pe++)etype(t,u,e[s]);
        e[t]=new(Pe++)etype(s,0,e[t]);
        e[s]->pair=e[t];
        e[t]->pair=e[s];
    }
    memmove(eb,e,sizeof(e));
    memset(h,0,sizeof(h));
    memset(vh,0,sizeof(vh));
    vh[0]=n;
    int ans=0;
    while(h[1]<n)ans+=aug(1,~0U>>1);
    printf("%d\n",ans);
    return 0;
}

This post has been edited 2 times. Last edited by zkw, Jul 6, 2011, 5:35 PM

SPFA 的两个优化

by zkw, Jul 6, 2011, 5:14 PM

　　SPFA 与堆优化的 Dijkstra 的速度之争不是一天两天了, 不过从这次 USACO 月赛题来看, SPFA 用在分层图上会比较慢. 标程是堆优化的 Dijkstra, 我写了一个非常朴素的 SPFA, 只能过 6/11 个点. SPFA 是按照 FIFO 的原则更新距离的, 没有考虑到距离标号的作用. 实现中 SPFA 有两个非常著名的优化: SLF 和 LLL.

　　SLF: Small Label First 策略.
　　实现方法是, 设队首元素为 $i$ , 队列中要加入节点 $j$ , 在 $d_j \le d_i$ 时加到队首而不是队尾, 否则和普通的 SPFA 一样加到队尾.

　　LLL: Large Label Last 策略.
　　实现方法是, 设队列 $Q$ 中的队首元素为 $i$ , 距离标号的平均值为 $\overline d = \frac {\sum_{j \in Q} d_j }{\left| Q \right|}$ , 每次出队时, 若 $d_i > \overline d$ , 把 $i$ 移到队列末尾, 如此反复, 直到找到一个 $i$ 使 $d_i \le \overline d$ , 将其出队.

　　加上 SLF 优化后程序多了一行, 过了 9/11 个点. 你问我怎么用 SPFA AC 这个题？利用分层图性质, 算完一层再算一层, 对每一层计算用 SPFA, 加上上面的优化, 程序飞快: 最强的优化要利用题目的特殊性质.

Quote:

大家怎么看?

This post has been edited 1 time. Last edited by zkw, Mar 22, 2014, 5:30 PM

从入门到精通: 最小费用流的“zkw算法”

by zkw, Jul 6, 2011, 4:48 PM

1. 网络流的一些基本概念 点击阅读

2. 最小费用流的各种转化 点击阅读

3. 费用流中的负边和负圈 点击阅读

　　在费用流的求解中难免会遇到负边和负圈的问题. 对于消圈算法, 负圈就是算法本身的一部分; 但对于增广路算法, 负圈是我们所不愿意看到的. 鉴于竞赛中使用消圈算法将带来严重的效率问题, 我们必须探索使用增广路算法处理费用流负边和负圈的方法.

　　对于单纯的负边, 如果这些负边没有组成负圈, 可以使用重赋权技术来处理, 即通过对每个节点合适的顶标, 使得 Reduced Cost 非负. 这个顶标通常可以使用到汇点的距离, 修改之后的边权变为 $c_{ij}^\pi = c_{ij} - D_i + D_j$ . 根据流量平衡条件, 容易根据这个新的费用算出原费用, 同时可以证明这样的重赋权不改变最优解. 这样做的代价是一次 SPFA 操作, 时间耗费是相对较小的.

　　如果存在负圈, SPFA 算法将不会终止. 很多同学使用人为的限制条件使其终止, 这是错误的. 举一个例子: 源点和汇点均为孤立点, 图中另外存在一个负费用, 正容量的圈. 不管怎样初始标号, 仅凭增广无法消除这个负圈, 从而得不到最优解 (根据最小费用最大流的定义, 要在所有最大流 (流量都为 $0$ ) 中寻找一个费用最小的方案). 也许你会说这个例子太过极端, 但是也很容易构造出其他一些例子, 说明不处理负圈, 或者简单地对算法打补丁并不能解决本质问题.

　　解决方法就是将所有负费用的边强制满流, 称为“退流”操作. 这样做之后我们破坏了平衡条件, 但满足了最优条件. 之后我们可以先求可行流, 再求最大流, 其间一直维持最优条件, 并逐步解决平衡条件, 最大流条件等问题. 这也就是 (2.) 中提到的方法. 这个方法可以消除所有负权边 (在 Reduced Cost 意义下), 同时正确处理所有负圈问题. 那么, 这个方法的速度如何呢?

　　费用流相关的图大致可以分为两类, 一类图侧重于费用, 即总的流量不大, 如 Two Shortest (只有 $2$ ), KaKa's Matrix Travels 等. 较为通用的强制满流方法在这类图上表现不佳, 原因是原本的流量较小而负费用较多, 经过转化的图在可行流求解阶段流量与原流量相比很大, 严重影响速度. 另一方面, 这类图往往有深刻的最短路背景, 从而不会出现负圈, 可以使用 SPFA 重赋权. 另一类图侧重于流量, 即边的费用相差不大, 如 employee, 最优图像的求解. 这种图中流量不小而负费用相对有限 (也可以稍大, 关键是前者), 经过转化增加的流量可以很快完成计算, 因此强制满流方法就没有问题.

　　另外, 不同的建图方式有可能造成不同的模型. 在 employee 这道题中, 如果建图时从每个点向后连边, 从前向后运行最小费用最大流, 这时的算法就只有负边, 没有负圈. 而如果从每个点向前连边, 在整个图上求最小费用平衡流, 就会有负圈出现. 但是这里并没有本质区别: 将第二种模型中的所有向前连的边强制满流即可得到第一种模型, 如果将其他的负边也强制满流还能得到一个完全没有负边的模型. 大家可能也猜到了, 这几种模型在速度上相去不远, 毕竟根据我们刚才的讨论, 增加的流量在算法的执行中不占主要地位.

4. 最小费用流的“zkw算法” 点击阅读

5. “zkw” 费用流算法在哪些图上慢 点击阅读

　　实践中, 上面的这个算法非常奇怪. 在某一些图上, 算法速度非常快, 另一些图上却比纯 SPFA 增广的算法慢. 不少同学经过实测总结的结果是稠密图上比较快, 稀疏图上比较慢, 但也不尽然. 这里我从理论上分析一下, 究竟这个算法用于哪些图可以得到理想的效果.

　　先分析算法的增广流程. 和 SPFA 直接算法相比, 由于同属于沿最短路增广的算法, 实际进行的增流操作并没有太多的区别, 每次的增流路径也大同小异. 因此不考虑多路增广时, 增广次数应该基本相同. 运行时间上主要的差异应当在于如何寻找增流路径的部分.

　　那么 zkw 算法的优势在于哪里呢? 与 SPFA 相比, KM 的重标号方式明显在速度上占优, 每次只是一个对边的扫描操作而已. 而 SPFA 需要维护较为复杂的标号和队列操作, 同时为了修正标号, 需要不止一次地访问某些节点, 速度会慢不少. 另外, 在 zkw 算法中, 增广是多路进行的, 同时可能在一次重标号后进行多次增广. 这个特点可以在许多路径都费用相同的时候派上用场, 进一步减少了重标号的时间耗费.

　　下面想一想 zkw 算法的劣势, 也就是 KM 重标号方式存在的问题. KM 重标号的主要问题就是, 不保证经过一次重标号之后能够存在增广路. 最差情况下, 一次只能在零权网络中增加一条边而已. 这时算法就会反复重标号, 反复尝试增广而次次不能增广, 陷入弄巧成拙的境地.

　　接下来要说什么, 大家可能已经猜到了. 对于最终流量较大, 而费用取值范围不大的图, 或者是增广路径比较短的图 (如二分图), zkw 算法都会比较快. 原因是充分发挥优势. 比如流多说明可以同一费用反复增广, 费用窄说明不用改太多距离标号就会有新增广路, 增广路径短可以显著改善最坏情况, 因为即使每次就只增加一条边也可以很快凑成最短路. 如果恰恰相反, 流量不大, 费用不小, 增广路还较长, 就不适合 zkw 算法了.

　　不适合怎么办? 继续向下看.

6. 最小费用流的原始对偶 (Primal-Dual) 算法 点击阅读

　　最小费用流的直接 SPFA 算法和前面说的 KM 重标号算法, 各自都有一些情况非常慢. 这里要写的就是一个所谓的“新算法” (其实非常经典), 融合两者优势.

　　先从 zkw 算法的本质开始. 细心的同学已经注意到了, 算法的主要过程就是反复交替进行最短路和最大流的计算. 具体地说, 每一次增广过程就是在零权网络中找一个最大流增广, 而重标号过程就是找到一组新的可行的距离标号. 既然如此, 最大流和最短路就完全可以用其他方法来替换. 这种交替进行最短路和最大流计算用来求得最小费用最大流的算法, 其实就是非常经典的原始对偶算法.

　　费用流的算法大致分为两种, 一种是经典的解法, 如消圈, 增广路, 原始对偶等等, 特点是步步为营, 维持可行性或者最优性其中之一, 再不断对另一方面作出改进. 另一种就比较现代一些, 典型的例子是松弛算法和网络单纯形, 由于放松了对求解过程中解的限制条件, 使得其速度远远超过经典解法, 同时也增加了编程难度和理解障碍. 下面要说的原始对偶算法, 速度自然不可能比松弛和网络单纯形快, 但应该是经典解法中的佼佼者了.

　　先考虑最短路算法的选择. 在流大, 费用小, 距离短的图上, KM 重标号就可以取得不错的效果, 但是 SPFA 在另外的图上就很有优势.

　　再考虑最大流算法的选择. 最朴素的 SPFA 暴力增广中根本就不维护距离标号, 这里隐含地使用了一个单路增广. 在 zkw 费用流的算法里, 这里使用的是多路增广. 是不是有必要用更高级的网络流算法呢? 经过实际测试, 除非在特别稠密的图上, 否则效果一般. 原因是可增广的流量并不是很大, 如果去维护 SAP 的顶标反而增加了不少开销, 预流推进亦然.

　　总结一下, 折中的选择就是使用 Small Label First 优化的 SPFA 来维护 zkw 算法中的距离标号, 保留多路增广. 注意, 除了第一次之外, SPFA 算法都工作在一个所有边均为正权的图上 (Reduced cost 意义下), 这是和不维护顶标直接 SPFA 的不同之处, 也是原来 zkw 算法的精神所在. 也正因为如此, 这里的 SPFA 甚至可以用 Dijkstra 替换.

　　这个特殊的原始对偶算法在稠密二分费用小的图上不敌原来的 zkw 算法, 但远远胜过暴力 SPFA. 在另外的图上, 对两者都是稳胜. 比暴力 SPFA 快原因是, 多路增广, 同时使用了 Reduced Cost 缩小了费用范围, 从而利于 SPFA 算法的工作 (需要的松弛次数减少), 而且使用 Reduced Cost 后不再有负边, 使 SLF 的优化落到了实处 (回忆: SLF 优化只有当所有边均为正的时候才能发挥出最佳效果), 甚至允许用 Dijkstra 来完成后面的工作. 比 zkw 算法快的原因是, 在流小费用大距离长的图上, 一次性把距离标号改对往往比反复调整更有效率.

　　最后是代码, 这次的例子是 POJ 3680 zkw Accepted 216K 360MS. 另外还有一个 POJ 3422, 是原来的 zkw 算法不擅长的一个例子, 这种方法只用 47MS, 代码大同小异, 就不放在这里了.

　　点击查看程序代码

#include <cstdio>
#include <cstring>
#include <deque>
#include <algorithm>
using namespace std;

const int V=440, E=V*2, maxint=0x3F3F3F3F;

struct etype
{
    int t, c, u;
    etype *next, *pair;
    etype() {}
    etype(int T, int C, int U, etype* N): t(T), c(C), u(U), next(N) {}
    void* operator new(unsigned, void* p){return p;}
} *e[V], Te[E+E], *Pe;

int S, T, n, piS, cost;
bool v[V];

void addedge(int s, int t, int c, int u)
{
    e[s] = new(Pe++) etype(t, +c, u, e[s]);
    e[t] = new(Pe++) etype(s, -c, 0, e[t]);
    e[s]->pair = e[t];
    e[t]->pair = e[s];
}

int aug(int no, int m)
{
    if (no == T) return cost += piS * m, m;
    v[no] = true;
    int l = m;
    for (etype *i = e[no]; i; i = i->next)
        if (i->u && !i->c && !v[i->t])
        {
            int d = aug(i->t, l < i->u ? l : i->u);
            i->u -= d, i->pair->u += d, l -= d;
            if (!l) return m;
        }
    return m - l;
}

bool modlabel()
{
    static int d[V]; memset(d, 0x3F, sizeof(d)); d[T] = 0;
    static deque<int> Q; Q.push_back(T);
    while(Q.size())
    {
        int dt, no = Q.front(); Q.pop_front();
        for(etype *i = e[no]; i; i = i->next)
            if(i->pair->u && (dt = d[no] - i->c) < d[i->t])
                (d[i->t] = dt) <= d[Q.size() ? Q.front() : 0]
                    ? Q.push_front(i->t) : Q.push_back(i->t);
    }
    for(int i = 0; i < n; ++i)
        for(etype *j = e[i]; j; j = j->next)
            j->c += d[j->t] - d[i];
    piS += d[S];
    return d[S] < maxint;
}

int ab[V], *pab[V], w[V];

struct lt
{
    bool operator()(int* p1,int* p2) {return *p1 < *p2;}
};

int main()
{
    int t;
    scanf("%d",&t);
    while(t--)
    {
        memset(e,0,sizeof(e));
        Pe = Te;
        static int m, k;
        scanf("%d %d", &m, &k);
        int abz = 0;
        for(int i = 0; i < m; ++i)
        {
            scanf("%d", pab[abz] = &ab[abz]), abz++;
            scanf("%d", pab[abz] = &ab[abz]), abz++;
            scanf("%d", &w[i]);
        }
        sort(&pab[0], &pab[abz], lt());
        int c=0xDEADBEEF; n=0;
        for(int i = 0; i < abz; ++i)
        {
            if(c != *pab[i]) c = *pab[i], ++n;
            *pab[i] = n;
        }
        ++n, S = 0, T = n++;
        for(int i = 0; i < T; ++i) addedge(i, i+1, 0, k);
        for(int i = 0; i < m; ++i) addedge(ab[i+i], ab[i+i+1], -w[i], 1);
        piS = cost = 0;
        while(modlabel())
            do memset(v, 0, sizeof(v));
            while(aug(S, maxint));
        printf("%d\n", -cost);
    }
    return 0;
}

This post has been edited 8 times. Last edited by zkw, Jun 2, 2020, 8:15 AM
Reason: Typo

坐井观天