[Dr.Lib]Note:解题报告 – 上网（DP）

题4 上网

【问题描述】

假设有n个人要上网，却只有1台电脑可以上网。上网的时间是从1 szw 至 T szw ，szw是sxc,zsx,wl自创的时间单位，至于 szw怎么换算成s,min或h，没有人清楚。依次给出每个人在某个时间段内上网的快乐程度C（必须这个人在整个时间段内都在上网，才能获得快乐程度C，否则，快乐程度是0），请你得到使总的快乐程度达到最大的方案。

【输入格式】

第1行2个整数 n和T，含义如题目所述；

接下来有n个这样的结构（每两个相邻的结构之间有一空行，且第1个结构和第一行间有一空行）：

第1行一个整数Mi，表示第i个人的时间段的个数；

接下来有Mi行，每行3个整数Xj,Yj,C，表示第i个人在[Xj，Yj]内上网的快乐程度为C，

因此有Xj-Yj-1=1，X1=1，Ymi=T，Xj<=Yj。

【输出格式】

仅输出一行，为总的最大的快乐程度。

【输入样例】

3 10

1 3 6

4 7 9

8 10 3

1 3 5

4 7 10

8 10 1

1 3 2

4 8 2

9 9 6

10 10 3

【输出样例】

【样例说明】

在[1，3]内，安排1上网，快乐程度为6；

在[4，7]内，安排2上网，快乐程度为10；

在[8，8]内，不安排；

在[9，9]内，安排3上网，快乐程度为6；

在[10，10]内，安排3上网，快乐程度为3；

这是使总的快乐程度达到最大的方案，对应的值是25。

【数据范围】

对于30%的数据，n<=4，所有的Mi<=5，T<=20；

对于60%的数据，n<=100，所有的Mi<=100，T<=2000；

对于100%的数据，n<=500，所有的Mi<=500，T<=500000，所有的0<C<=10^9，并保证最终解Max<=10^9。

【标程】

很明显,前30%的数据可以用暴力过，直接枚举每一个区间即可，时间复杂度是2^(n*m)。

其实，拿到这到题，就很容易想到DP，用max[i]表示第i szw后，所能达到的最大快乐程度，状态转移方程如下：max[i]=max{max[i-1]，max[x-1]+C[x,i]}，边界情况max[0]=0，其中C[x,i]表示从x时刻开始玩到i时刻的快乐程度。具体操作时，枚举所有结束时刻是i的区间，比较这些区间的初始时刻x对应的max[x-1]+C[x,i]与当前max[i]的大小，如果max[x-1]+C[x,i]更大，就替换max[i]。可能对于很多选手来说，拿出状态转移方程并不难，问题在于怎么存每个区间。如果直接开2维数组，c:array[1..maxT,1..maxT] of longint，空间复杂度是T^2，可以过60%的数据。我们应该如何存区间呢？

我们开一个一维数组，g:array[1..maxn*maxm] of record x,y,c,l:longint; end，空间约4M，其中x,y表示对应区间的初始时刻和结束时刻，c表示快乐程度，l表示上一个和这个区间y值相同的区间所对应的位置。再开一个一维数组，l:array[1..maxT] of longint，l[y]表示以y为结束时刻的最后一个区间所对应的位置。这样，我们就可以很轻松的找到以某个值y为结束时刻的区间了。DP的时间复杂度降到 O(maxn*maxm+T)。但事实上由于读入大数据文件比较耗时,所以标程过最后一个数据也要约0.4秒，因此把时限定为1.5s。(我们同样可以用x做为标准)

此题主要考察动态规划和链表，拿60分相对容易，拿满分比较困难。

数据范围和类型:

数据	1	2	3	4	5	6	7	8	9	10
N值	3	4	4	65	78	98	400	450	498	500
M值	4	4	5	62	92	98	384	462	498	500
T值	10	16	19	1234	1674	98	25*10^4	35*10^4	45*10^4	5*10^5
平均分	3.81	2.90	3.24	2.33	2.33	2.84	1.19	1.19	1.14	1.08
类型	手工	随机	特殊	随机	随机	特殊	随机	随机	随机	随机

program hh;
var f:array[0..500010] of longint;
    first:array[1..500000] of longint;
    e:array[1..250000] of record x,y,next,v:longint; end;
    i,j,k,n,t,s,x,y,v:Longint;
function max(x,y:Longint):Longint;
begin
 if x>y then exit(x) else exit(y);
end;

begin
 assign(input,'net.in');
 reset(input);
 assign(output,'net.out');
 rewrite(output);
 readln(n,t);
 s:=0;
 for i:=1 to n do
  begin
   readln;
   readln(k);
   for j:=1 to k do
    begin
     inc(s);
     readln(x,y,v);
     e[s].x:=x;e[s].y:=y;e[s].v:=v;e[s].next:=first[y];first[y]:=s;
    end;
   end;
 for i:=1 to t do
  begin
   f[i]:=f[i-1]; k:=first[i];
   while k<>0 do
    begin
     f[i]:=max(f[e[k].x-1]+e[k].v,f[i]);
     k:=e[k].next;
    end;
  end;
 writeln(f[t]);
 close(input);close(output);
end.

program hh;

var f:array[0..500010] of longint;

first:array[1..500000] of longint;

e:array[1..250000] of record x,y,next,v:longint; end;

i,j,k,n,t,s,x,y,v:Longint;

function max(x,y:Longint):Longint;

begin

if x>y then exit(x) else exit(y);

end;

begin

assign(input,'net.in');

reset(input);

assign(output,'net.out');

rewrite(output);

readln(n,t);

s:=0;

for i:=1 to n do

begin

readln;

readln(k);

for j:=1 to k do

begin

inc(s);

readln(x,y,v);

e[s].x:=x;e[s].y:=y;e[s].v:=v;e[s].next:=first[y];first[y]:=s;

end;

for i:=1 to t do

begin

f[i]:=f[i-1]; k:=first[i];

while k<>0 do

begin

f[i]:=max(f[e[k].x-1]+e[k].v,f[i]);

k:=e[k].next;

end;

writeln(f[t]);

close(input);close(output);

end.

【个人报告】

赛场上前几题耗时太多没来得及打这题。妥妥的DP阿……

但确实没想那么多……

看了题解才知道确实有坑。

抽象的模型大概就是m条带权线段求不相交区间最大权值和……几个几个人是考你读入的……

果断读进来按终点排序，一条一条来DP。

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<vector>
#include<cmath>
#include<queue>
#include<stack>
using namespace std;
typedef long long ll;
int n,t,mt,m,px,py,pc;
struct net{
	int x;int y;int c;
};
struct p{
	int y;int c;
};
struct dqueue{
	dqueue(){head=-1;}
	p S[250000];
	int head;
	bool ins(p i){
		if(head==-1){head++;S[head]=i;return 1;}
		if(S[head].c>=i.c)return 0;
		head++;S[head]=i;return 1;
	}
	int get(int i){
		if(head==-1){return 0;}
		int t=head;
		while(t>=0&&S[t].y>i){--t;}
		if(t==-1)return 0;
		return S[t].c;//二分查找我不会= =
	}
};
int cmp(net a,net b){
	return a.y<b.y;
}
vector<net> Q;
dqueue D;
int scanfd(){
	int p=0;
	char c;
	while(c=getchar()){
		if(c-'0'<=9&&c-'0'>=0){break;}
	}
	p=c-'0';
	while(c=getchar()){
		if(c-'0'<=9&&c-'0'>=0){p=p*10+c-'0';}else{break;}
	}
	return p;
}
int main(){
    freopen("net.in","r",stdin);
    freopen("net.out","w",stdout);
	scanf("%d%d",&n,&t);
	while(n--){
		scanf("%d",&m);
		mt+=m;
		while(m--){
			px=scanfd();py=scanfd();pc=scanfd();
			net pn;
			pn.x=px;pn.y=py;pn.c=pc;
			Q.push_back(pn);
		}
	}
	sort(Q.begin(),Q.end(),cmp);
	for(int i=0;i<mt;++i){
		pc=Q[i].c;px=Q[i].x;py=Q[i].y;
		int tmc=D.get(px-1);
		p pyc;
		pyc.c=tmc+pc;pyc.y=py;
		D.ins(pyc);
	}
	cout<<D.get(t)<<endl;
	return 0;
}

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cstring>

#include<algorithm>

#include<vector>

#include<cmath>

#include<queue>

#include<stack>

using namespace std;

typedef long long ll;

int n,t,mt,m,px,py,pc;

struct net{

int x;int y;int c;

};

struct p{

int y;int c;

};

struct dqueue{

dqueue(){head=-1;}

p S[250000];

int head;

bool ins(p i){

if(head==-1){head++;S[head]=i;return 1;}

if(S[head].c>=i.c)return 0;

head++;S[head]=i;return 1;

}

int get(int i){

if(head==-1){return 0;}

int t=head;

while(t>=0&&S[t].y>i){--t;}

if(t==-1)return 0;

return S[t].c;//二分查找我不会= =

}

};

int cmp(net a,net b){

return a.y<b.y;

}

vector<net> Q;

dqueue D;

int scanfd(){

int p=0;

char c;

while(c=getchar()){

if(c-'0'<=9&&c-'0'>=0){break;}

}

p=c-'0';

while(c=getchar()){

if(c-'0'<=9&&c-'0'>=0){p=p*10+c-'0';}else{break;}

}

return p;

}

int main(){

freopen("net.in","r",stdin);

freopen("net.out","w",stdout);

scanf("%d%d",&n,&t);

while(n--){

scanf("%d",&m);

mt+=m;

while(m--){

px=scanfd();py=scanfd();pc=scanfd();

net pn;

pn.x=px;pn.y=py;pn.c=pc;

Q.push_back(pn);

}

sort(Q.begin(),Q.end(),cmp);

for(int i=0;i<mt;++i){

pc=Q[i].c;px=Q[i].x;py=Q[i].y;

int tmc=D.get(px-1);

p pyc;

pyc.c=tmc+pc;pyc.y=py;

D.ins(pyc);

}

cout<<D.get(t)<<endl;

return 0;

}

[Dr.Lib]Note:Math – 不定方程

Via Ivy-Endhttp://www.ivy-end.com/archives/1010
关于这个算法，主要是参考NOIP2012 Day2 T1。即这里所讲的是求解这样一个线性模方程：\[ax\equiv 1\mod{p}\]的最小正整数解。

去年我是暴搜做的(PS+1)，当时什么都不会。今天在这里介绍两种算法，一种是我国古代数学家秦九韶发明的「大衍求一术」，还一种是著名的「扩展欧几里德算法」。

大衍求一术

首先来看一下大衍求一术。这里只介绍它的计算方法，至于证明可以参考扩展欧几里德算法。

例1：求解方程\(23x\equiv 1\mod{97}\)。

解：我们只需要列出下面这张表就可以得到求解\[\begin{matrix}23^{1} & 23^{1} & 3^{17} & 3^{17} & 1^{38}\\ 97^{0} & 5^{4} & 5^{4} & 2^{21} & 2^{21}\end{matrix}\]结果就是38。

接下来我们来理论化的表述一下这个算法的过程：

假设输入\(a,b\)满足\(a>b\)。那么我们用\(a_{n},A_{n}\)分别表示第一行的底数和奇数，\(b_{n},B_{n}\)分别表示第二行的底数和奇数，如果\(a_{i}>b_{i}\)，那么\(a_{i+1}=a_{i}\mod{b_{i}},A_{i+1}=A_{i}+B_{i}\cdot \left [ \frac{a_{i}}{b_{i}} \right ],b_{i+1}=b_{i},B_{i+1}=B_{i}\)；如果\(a_{i}<b_{i}\)则上面的结论倒过来即可。

算法结束当且仅当\(a_{i}=1\)，此时\(A_{i}\)即为所求的最小正整数解。

例2：求解方程\(97x\equiv 1\mod{23}\)。

解：我们只需要列出下面这张表就可以得到求解\[\begin{matrix}97^{1} & 5^{1} & 5^{1} & 2^{5} & 2^{5} & \\ 23^{0} & 23^{0} & 3^{4} & 3^{4} & 1^{9} & 1^{14}\end{matrix}\]结果就是14。

对于这个结果，如果1最先出现在下面一行，则需要再计算一次，而且这次计算必须使得余数是1。

假设输入\(a,b\)满足\(a<b\)。中间的步骤和之前一行，在计算过程中必然存在一个\(i\)使得\(b_{i}=1\)，此时我们只需计算\(B_{i+1}\)即可得到结果。其中\(B_{i+1}=A_{i}+B_{i}\cdot \left(a_{i} – 1\right)\)。

#include <iostream>

using namespace std;

struct Num
{
    int nBase, nIndex;    // 分别表示底数和奇数
};

Num x, y;

int main()
{
    cin >> x.nBase >> y.nBase;
    x.nIndex = 1; y.nIndex = 0;    // 初始化
    if(x.nBase < y.nBase)
    {
        while(1)
        {
            if(x.nBase == 1 || x.nBase == 0)    // 循环出口
            { cout << x.nIndex << endl; break; }
            if(x.nBase < y.nBase)    // 模拟计算过程
            {
                int nDiv = y.nBase / x.nBase;
                if(y.nBase % x.nBase == 0)
                { nDiv--; }
                y.nBase %= x.nBase;
                y.nIndex += nDiv * x.nIndex;
            }
            else
            {
                int nDiv = x.nBase / y.nBase;
                if(x.nBase % y.nBase == 0)
                { nDiv--; }
                x.nBase %= y.nBase;
                x.nIndex += nDiv * y.nIndex;
            }
        }
    }
    else
    {
        while(1)
        {
            if(y.nBase == 1)    // 出口
            {
                int nDiv = x.nBase - 1;
                cout << nDiv * y.nIndex + x.nIndex << endl;    // 求出结果
                break;
            }
            if(x.nBase < y.nBase)    // 模拟计算过程
            {
                int nDiv = y.nBase / x.nBase;
                if(y.nBase % x.nBase == 0)
                { nDiv--; }
                y.nBase %= x.nBase;
                y.nIndex += nDiv * x.nIndex;
            }
            else
            {
                int nDiv = x.nBase / y.nBase;
                if(x.nBase % y.nBase == 0)
                { nDiv--; }
                x.nBase %= y.nBase;
                x.nIndex += nDiv * y.nIndex;
            }
        }
    }
    return 0;
}

#include <iostream>

using namespace std;

struct Num

{

int nBase, nIndex; // 分别表示底数和奇数

};

Num x, y;

int main()

{

cin >> x.nBase >> y.nBase;

x.nIndex = 1; y.nIndex = 0; // 初始化

if(x.nBase < y.nBase)

{

while(1)

{

if(x.nBase == 1 || x.nBase == 0) // 循环出口

{ cout << x.nIndex << endl; break; }

if(x.nBase < y.nBase) // 模拟计算过程

{

int nDiv = y.nBase / x.nBase;

if(y.nBase % x.nBase == 0)

{ nDiv--; }

y.nBase %= x.nBase;

y.nIndex += nDiv * x.nIndex;

}

else

{

int nDiv = x.nBase / y.nBase;

if(x.nBase % y.nBase == 0)

{ nDiv--; }

x.nBase %= y.nBase;

x.nIndex += nDiv * y.nIndex;

}

else

{

while(1)

{

if(y.nBase == 1) // 出口

{

int nDiv = x.nBase - 1;

cout << nDiv * y.nIndex + x.nIndex << endl; // 求出结果

break;

}

if(x.nBase < y.nBase) // 模拟计算过程

{

int nDiv = y.nBase / x.nBase;

if(y.nBase % x.nBase == 0)

{ nDiv--; }

y.nBase %= x.nBase;

y.nIndex += nDiv * x.nIndex;

}

else

{

int nDiv = x.nBase / y.nBase;

if(x.nBase % y.nBase == 0)

{ nDiv--; }

x.nBase %= y.nBase;

x.nIndex += nDiv * y.nIndex;

}

return 0;

}

扩展欧几里德算法

可能上面的算法对于某些人来说比较晦涩，我们下面来介绍一下扩展欧几里德算法。首先介绍一个定理：

方程\(ax+by=\gcd\left ( a,b \right )\)一定有解。

这样我们的问题就可以转化为求方程\(ax+b\cdot \left ( -y \right )=1\)，在这里，我们先求出方程\(ax+b\cdot \left ( -y \right )=\gcd\left(a,b\right)\)的解，然后只要将结果除以\(\gcd\left(a,b\right)\)就行了。

下面来推导一下扩展欧几里德算法。

我们已知\[ax+by=\gcd\left ( a,b \right )\]，且\(\gcd\left ( a,b \right )=\gcd\left(b,a\mod b \right )\)。不妨设\[bx{}'+\left ( a\mod b \right )y{}'=\gcd\left ( b,a\mod b \right )\]。此时就有\[bx{}'+\left ( a\mod b \right )y{}'=ax+by\]，展开得到\[bx{}'+\left ( a-\left [ \frac{a}{b} \right ]\cdot b \right )y{}’=ax+by\]，化简得\[ay{}'+b\left (x{}'-\left [ \frac{a}{b} \right ]\cdot y{}’ \right )=ax+by\]。因此可以得到\[x=y{}',y=x{}'-\left [ \frac{a}{b} \right ]\cdot y{}’\]。

这样我们就可以用递归来实现扩展欧几里德算法了。

#include <iostream>

using namespace std;

typedef long long LL;

LL A, B, C, X = 0, Y = 0;

LL gcd(LL a, LL b);
void exgcd(LL a, LL b, LL &x, LL &y);

int main()
{
    cin >> A >> B;
    C = gcd(A, B);
    exgcd(A, -B, X, Y);
    while(X < 0) { X += B; }    // 找最小正整数
    cout << X << endl;
    return 0;
}

void exgcd(LL a, LL b, LL &x, LL &y)
{
    if(a == 0)
    {
        x = 0; y = C / b;    // 边界
    }
    else
    {
        exgcd(b % a, a, x, y);
        y = x;        // 递推公式
        x = (C - b * y) / a;
    }
}

LL gcd(LL a, LL b)    // 求解最大公倍数
{
    if(b == 0) { return a; }
    else
    {
        return gcd(b, a % b);
    }
}

#include <iostream>

using namespace std;

typedef long long LL;

LL A, B, C, X = 0, Y = 0;

LL gcd(LL a, LL b);

void exgcd(LL a, LL b, LL &x, LL &y);

int main()

{

cin >> A >> B;

C = gcd(A, B);

exgcd(A, -B, X, Y);

while(X < 0) { X += B; } // 找最小正整数

cout << X << endl;

return 0;

}

void exgcd(LL a, LL b, LL &x, LL &y)

{

if(a == 0)

{

x = 0; y = C / b; // 边界

}

else

{

exgcd(b % a, a, x, y);

y = x; // 递推公式

x = (C - b * y) / a;

}

LL gcd(LL a, LL b) // 求解最大公倍数

{

if(b == 0) { return a; }

else

{

return gcd(b, a % b);

}

欧拉定理

\(a^{phi(n)}=1 \mod {n}\)
令\(x=a^{phi(n) -1}\mod{n}\)有\(ax=1 \mod{n}\)则x即为答案。
只需求出phi(n)就可以了。

#include<iostream>
#include<cmath>
using namespace std;
typedef unsigned int ll;
typedef long long lll;
lll qkpow(lll a,  lll b,  lll  c){
    lll  ret=1%c;
   while(b) {
      if(b&1) ret=ret*a%c;a=a*a%c;b>>=1;
   }
   return ret;
}
ll a,b,n;lll rsts;
ll phi(){
    ll rst=1;
    double x=sqrt(b);
    for(int i=2;i<=x;++i){
        if(b%i==0){
            rst*=i-1;
            while(b%i==0){
                rst*=i;
                b/=i;
            }
            rst/=i;
        }
    }
    if(b>1){rst*=(b-1);}
    return rst;
}
int main(){
    cin>>a>>b;n=b;
    rsts=qkpow(a,phi()-1,n);
    cout<<rsts<<endl;
    return 0;
}

#include<iostream>

#include<cmath>

using namespace std;

typedef unsigned int ll;

typedef long long lll;

lll qkpow(lll a, lll b, lll c){

lll ret=1%c;

while(b) {

if(b&1) ret=ret*a%c;a=a*a%c;b>>=1;

}

return ret;

}

ll a,b,n;lll rsts;

ll phi(){

ll rst=1;

double x=sqrt(b);

for(int i=2;i<=x;++i){

if(b%i==0){

rst*=i-1;

while(b%i==0){

rst*=i;

b/=i;

}

rst/=i;

}

if(b>1){rst*=(b-1);}

return rst;

}

int main(){

cin>>a>>b;n=b;

rsts=qkpow(a,phi()-1,n);

cout<<rsts<<endl;

return 0;

}

[Dr.Lib]Note:Algorithms – 矩阵乘法

矩阵乘法在很多题目中都有显著的应用，在一些递推、计算几何和其他可以用矩阵表示元素操作的题目中，可以借助快速幂的思想优化不少时间。

1.加速递推

最著名的递推当属Fibonacci数列了。
\(
f[0]=f[1]=1
f[n]=f[n-1]+f[n-2]
\)

朴素的递推需要O(n)的时间（应该没人写递归吧……

如果采用矩阵快速幂

\(
\begin{bmatrix}
f(n+1) \\
f(n)
\end{bmatrix}
=
(\begin{bmatrix}
1 & 1\\
1 & 0
\end{bmatrix})^{n}
\cdot
\begin{bmatrix}
1 \\
1
\end{bmatrix}
\)
其中矩阵乘法可以采用快速幂的思想加速运算（矩阵乘法满足结合律）

稍稍复杂的例子还有

http://wikioi.com/problem/1281/

http://acm.buaa.edu.cn/problem/653/

……等等很多

2.计算几何

example：via http://www.matrix67.com/blog/archives/276/

给定n个点，m个操作，构造O(m+n)的算法输出m个操作后各点的位置。操作有平移、缩放、翻转和旋转
这里的操作是对所有点同时进行的。其中翻转是以坐标轴为对称轴进行翻转（两种情况），旋转则以原点为中心。如果对每个点分别进行模拟，那么m个操作总共耗时O(mn)。利用矩阵乘法可以在O(m)的时间里把所有操作合并为一个矩阵，然后每个点与该矩阵相乘即可直接得出最终该点的位置，总共耗时O(m+n)。假设初始时某个点的坐标为x和y，下面5个矩阵可以分别对其进行平移、旋转、翻转和旋转操作。预先把所有m个操作所对应的矩阵全部乘起来，再乘以(x,y,1)，即可一步得出最终点的位置。

3.……以及其他

example：黑书p205细菌、VOJ1049序列置换、有向图方案数……

顺次给出m个置换，反复使用这m个置换对初始序列进行操作，问k次置换后的序列。m<=10, k<2^31。
首先将这m个置换“合并”起来（算出这m个置换的乘积），然后接下来我们需要执行这个置换k/m次（取整，若有余数则剩下几步模拟即可）。注意任意一个置换都可以表示成矩阵的形式。例如，将1 2 3 4置换为3 1 2 4，相当于下面的矩阵乘法：

置换k/m次就相当于在前面乘以k/m个这样的矩阵。我们可以二分计算出该矩阵的k/m次方，再乘以初始序列即可。做出来了别忙着高兴，得意之时就是你灭亡之日，别忘了最后可能还有几个置换需要模拟。

给定一个有向图，问从A点恰好走k步（允许重复经过边）到达B点的方案数mod p的值
把给定的图转为邻接矩阵，即A(i,j)=1当且仅当存在一条边i->j。令C=A*A，那么C(i,j)=ΣA(i,k)*A(k,j)，实际上就等于从点i到点j恰好经过2条边的路径数（枚举k为中转点）。类似地，C*A的第i行第j列就表示从i到j经过3条边的路径数。同理，如果要求经过k步的路径数，我们只需要二分求出A^k即可。

4.模版

模版代码来源于http://blog.csdn.net/u010153200/article/details/8974369

template <int A, int B>
struct matrix
{
	long long v[A][B];
	matrix(long long r = 0)
	{
		memset(v, 0, sizeof(v));
		for (int i = 0; i < A && i < B; i++)
			v[i][i]	= r;
	}
	matrix(long long a[A][B])
	{
		memcpy(v, a, sizeof(a));
	}	
};
inline long long quickadd(long long a, long long b)
{
	long long ans = 0;
	for (; b; b >>= 1)
	{
		if (b & 1) ans = (ans + a) % m;
		a = (a + a) % m;
	}
	return ans;
}
template <int A, int B, int C>
matrix <A, C> operator * (const matrix <A, B> &a, const matrix <B, C> &b)
{
	matrix <A, C> ans;
	for (int i = 0; i < A; i++)
		for (int j = 0; j < C; j++)
		{
			long long v = 0;
			for (int k = 0; k < B; k++)
				v += quickadd(a.v[i][k], b.v[k][j]);
			ans.v[i][j] = (v % m + m) % m;
		}
	return ans;
}
template <int A>
matrix <A, A> power(const matrix<A, A> &a, long long b)
{
	matrix <A, A> ans = 1;
	for (matrix<A, A> i = a; b; b >>= 1)
	{
		if (b & 1) ans = ans * i;
		i = i * i;
	}
	return ans;
}

template <int A, int B>

struct matrix

{

long long v[A][B];

matrix(long long r = 0)

{

memset(v, 0, sizeof(v));

for (int i = 0; i < A && i < B; i++)

v[i][i] = r;

}

matrix(long long a[A][B])

{

memcpy(v, a, sizeof(a));

}

};

inline long long quickadd(long long a, long long b)

{

long long ans = 0;

for (; b; b >>= 1)

{

if (b & 1) ans = (ans + a) % m;

a = (a + a) % m;

}

return ans;

}

template <int A, int B, int C>

matrix <A, C> operator * (const matrix <A, B> &a, const matrix <B, C> &b)

{

matrix <A, C> ans;

for (int i = 0; i < A; i++)

for (int j = 0; j < C; j++)

{

long long v = 0;

for (int k = 0; k < B; k++)

v += quickadd(a.v[i][k], b.v[k][j]);

ans.v[i][j] = (v % m + m) % m;

}

return ans;

}

template <int A>

matrix <A, A> power(const matrix<A, A> &a, long long b)

{

matrix <A, A> ans = 1;

for (matrix<A, A> i = a; b; b >>= 1)

{

if (b & 1) ans = ans * i;

i = i * i;

}

return ans;

}

[Dr.Lib]Note:Algorithms – LLRB——红黑树的现代实现

Via http://www.cnblogs.com/Seiyagoo/archive/2013/10/13/3366290.html

一、本文内容

以一种简明易懂的方式介绍红黑树背后的逻辑实现2-3-4树，以及红黑树的插入、删除操作，重点在2-3-4树与红黑树的对应关系上，并理清红黑树相关操作的来龙去脉。抛弃以往复杂的实现，而分析红黑树的一种简单实现LLRB。

二、算法应用

红黑树，给人以强烈的第一听觉冲击力——红与黑，好像很高端的感觉。事实上的确如此，红黑树是一种高级数据结构，在C++、Java的标准库里作为set、map的底层数据结构实现，以及linux中进程的公平调度。

三、2-3-4树

标题是红黑树，为什么讲2-3-4树？因为红黑树就是2-3-4树的一种等价形式，更准确地来说，我们用红黑树来完成2-3-4树的各种操作（如插入、删除）。原因就是2-3-4树的实现即维护太麻烦。所以理解2-3-4树才能真正理解红黑树。而历史就是这么发展的，了解过去，现在的一切才有了意义。算法导论关于红黑树这一节就忽略了这一点，让人知其然而不知其所以然。

OK，暂时先忽略复杂的红黑树，从简单的2-3-4树开始。

1、定义

2-3-4树是一种泛化的BST，它的每个结点允许1,2或者3个键(key)，那么对应的有三种结点：

2-node：一个key，两个孩子；

3-node：二个key，三个孩子；

4-node：三个key，四个孩子。

注：k-node表示有k个链接（link）。泛化的BST还有2-3树，B树等。

从图中可以看出2-3-4树的另一个性质：它是完全平衡的（等高），即从根结点到叶子结点距离相等。

2、插入操作

2-3-4树本身就是一种查找树（中序遍历有序），故其查找操作同二叉查找。

2-3-4树的插入操作类似二叉查找树，先是查找操作失败（从根结点查找到叶子结点），然后在底部的叶子结点插入。

因为2-3-4树的结点有三种类型，所以操作有点差异。对于2-node和3-node，分别直接插入可变成3-node，4-node；但是对于4-node若直接插入则违反了定义。在4-node插入之前，先分裂4-node成2个2-node，再将待插入的key插入对应的2-node。如下图，H查找失败，在H插入4-node（由三个key F、G、J组成）之前，先对该4-node分裂（将三个key的中间值提上父节点，剩余的二个key分别作为中间key的左右孩子），然后再将H插入2-node J中。这样操作的结果是查找到达底部叶子结点时，始终是2-node或者3-node。

插入算法思想：自下而上的算法由原作者Bayer在1972年提出，自上而下的算法由Guibas-Sedgewick（红黑树这个名字来源于他们）在1978年提出，然后30年后也就是2008年Sedgewick教授又改进了红黑树的操作，也就是后面要介绍的LLRB。

自上而下的算法思路是，从根结点向下的查找过程中，遇到4-node就分裂，最后在底部的叶子结点插入。

那么为什么遇到4-node就分裂呢？4-node不是2-3-4树的一种合法结点类型吗？

答案可以从后面LLRB的算法思路可以得出。

因为遇到4-node就分裂就保证了当前结点不是4-node，则分裂孩子的4-node有两种情形：

分裂4-node的case 1

分裂4-node的case 2

注：上面的变换在树中任意位置都成立。

下面两张图是完整的插入过程（只有分裂结点类型为4-node的根结点才会导致树高增1）：

3、平衡性分析

2-3-4树的树高在最坏情况下为lgN（所有结点都是2-node型），最好情况下为lg4 N = 1/2 lgN（所有结点都是4-node型），2-3-4树的查找、插入操作都是lgN。

四、红黑树

终于到了高富帅——红黑树。。。

从2-3-4树的介绍可以看出，对2-node、3-node、4-node的不同数据类型进行转换，但所涉及的大部分任务使用这种直接的表示方法来实现并不方便。所以可以用

一种统一的方式完成转换，而只需很小的开销。这就是红黑树存在的意义，既有BST的标准搜索过程，又有2-3-4树的简单插入平衡过程。

下面介绍LLRB（Left-leaning red-black trees），而不是标准的红黑树。

1、定义

LLRB有三个特点：

（1）用BST来表示2-3-4树；

（2）用红边(红链接)来连接2-node来表示3-node和4-node（如下图）；

（3）3-node必须是向左倾斜的（两者的大者作为根）。

LLRB相对于标准的RB多了特点3，在标准的RB中右向倾斜的红链接是允许的。对于特点2，在物理上用一个bit（红或黑）来存储以表示指向该结点的红链接。

红链接来连接3-node或者4-node的内部key，而黑链接则连接外部的key；为了理解，可以消除红链接并将它们连接的结点都折叠起来（即将看做红链接连

接的点缩为一个点），则可以看出黑链接个数不变。

2-3-4树与红黑树是一一对应的关系

且上下关系中不允许2个连续的红边

由特点3可以推出LLRB的一个特性，红黑树与2-3-4树一一对应。

2、插入算法

同样地，在LLRB中查找操作同BST。

在插入之前要知道一个操作：旋转。它有两种情况：左旋，右旋。

左旋右旋

插入算法思路：即前面介绍的2-3-4树

具体实现时，插入一个结点时，始终是红结点，即用红边链接该结点。对于2-node、3-node直接插入（k-node有k个插入点），如违反上面的左红链接和连续的红链接，则旋转作调整。对于4-node（左右都为红链接），先分裂，物理实现是一个翻转（左右红链接变黑，父链接变红）。

2-node插入的两种case

3-node插入的三种case

4-node分裂操作

由4-node的分裂可知黑高度不变，分裂操作即翻转在图片上对应为红链接向上传递。

在介绍2-3-4树时，4-node分裂操作有两种情况，4-node的parent是2-node和3-node；再结合k-node有k个插入点，则总共有6种情况。

4-node的分裂case 1

4-node的分裂case 2

看了上面两幅图后，也许会让人觉得红黑树太复杂了，这么多case，其实不然，在LLRB实现中只有两种操作：旋转、翻转。旋转的目的是保持平衡，翻转的目的是分裂4-node。

看了下面的LLRB插入算法，你就会明白上面4-node的翻转、旋转其实是分开的两个过程（翻转自上而下，旋转自下而上），只是为了统一这个完整的过程而画在了一起，才会有那么多case。

LLRB的插入算法：

首先结合2-3-4树的插入算法思路，先从上至下查找（遇到4-node则翻转），然后在底部叶子结点插入，因为在从上至下的过程中，可能会产生不满足LLRB的性质的情况，故插入结点后需要从下至上调整以恢复LLRB性质。

下图是插入算法的核心代码，第2是分裂即翻转，第1是插入操作，第3、4是调整。

从插入算法可以看出，如果自下而上再分裂4-node，则会出现它的parent也可能是4-node，祖父结点也可能是4-node；我们可以一直向上分裂，这也正是上面提到的自上而下的思路；而更简单的方法是，在沿树向下的过程中，遇到4-node就分裂，这也正是自上而下与自下而上的区别。

插入算法的核心代码

上图的核心代码按照从上而下和从下而上的顺序放入BST的插入（递归版本）操作中即得到下图的完整的插入算法。

注：分裂（即翻转）是自上而下，所以放在递归之前；调整（即旋转）是自上而下，所以放在递归之后。

完整的插入代码

如果将分裂操作放到递归之后，也就是先自上而下查找，插入结点，然后自下而上调整也可同样完成插入操作而不破坏LLRB的性质。

2-3树的插入操作

其实上述描述的就是2-3树的插入操作，它与2-3-4树的插入的区别在于：2-3树先插入，再分裂（down）、调整（up）；2-3-4树先分裂（down），再插入、调整（up）。又因为插入总是在最后一层进行，故翻转的位置决定了对应树的实现。

这也是为什么2-3-4树叫top-down，而2-3树叫bottom-up。

3、删除算法

LLRB的删除类似于插入，只不过处理刚好相反。插入、删除都有临界点：插入4-node，删除2-node，对临界点的操作都会引起突变，因为它们会破坏LLRB的性质(黑高度)。

所以同插入一样，先从上至下查找，如果查找在3-node或4-node结束，则直接删除；

3-node和4-node的删除

对于2-node的删除同4-node的插入相反，2-node的删除是先合并2个2-node为1个4-node，然后再安全地删除对应的2-node中的key。

同样地，因为parent不为2-node（遇到即合并），再结合兄弟结点的2、3、4-node，则删除总共有6种情况。同样，实际的删除实现也

没这么复杂。

2-node的删除

在介绍删除任意一个结点时，先分析删除树中最小的结点。因为它是删除任意结点的一部分，后面可以看出来。

首先，为了保证可以直接删除最小的某个结点，需要假设当前结点h或者h.left是红色链。

然后从上而下查找过程中，2个2-node要变为1个4-node，则需反向翻转（红色父链接变黑，黑色子链接变红），

为了将红链从上向左子树传递（删除红结点，不改变黑高度），需保证h为红，h.left和h.left.left为黑；

当h.left和h.left.left都为黑时，

如果h.right.left为红，则要从右边借兄弟（下图case 2），如果h.right.left为黑，则不需要（下图case1）。

注：在翻转的同时，右子树可能会产生连续的红链，则需调整。

case 1

case 2

红链向左移动红链向左移动对应的example

deleteMin的实现

deleteMin的example

完成了deleteMin就完成了LLRB的删除操作的一大半。现在是删除LLRB的任意一个key，

自上而下查找过程中，左边查找用moveRedLeft；右边查找用moveRedRight；直到最后的底部叶子结点，直接删除即可；同样，自下而上调整。

怎样将delete操作归约到delteMin去呢？算法导论提供的一个技巧是：replace，deleteMin（即用后继的key代替当前的key，再删除右孩子的最小结点）。

删除技巧

完整删除代码

参考：

《算法导论》

《algorithm in c》

http://zh.wikipedia.org/wiki/%E7%BA%A2%E9%BB%91%E6%A0%91

http://www.cs.princeton.edu/~rs/talks/LLRB/RedBlack.pdf

[Dr.Lib]Note:Algorithms – 四边形不等式

若w[a,c]+w[b,d] \(\leq\) w[b,c]+w[a,d] (a<b<c<d)

则w满足凸四边形不等式

若w[i,j] \(\leq\) w[i,j] ([i,j]\(\subseteq\) [i',j'])

则w满足区间单调

2D/1D:d[i,j]=opt{d[i,k-1]+d[k,j]}+w[i,j] (1 \(\leq\) i < k \(\leq\) j \(\leq\) n)

d满足凸四边形不等式 k[i,j-1] \(\leq\) k[i,j] \(\leq\) k[i+1,j]

w为凸 \(\Leftrightarrow\) w[i,j]+w[i+1,j+1] \(\leq\) w[i+1,j]+w[i,j+1]

1D/1D:f[i]=opt{f[j]+w[j,i]} (b[i] \(\leq\) j \(\leq\) i-1)

满足决策单调性，可用单调队列优化决策

Tag Archives: algorithm

[Dr.Lib]Note:解题报告 – 上网（DP）

[Dr.Lib]Note:Math – 不定方程

大衍求一术

扩展欧几里德算法

欧拉定理

[Dr.Lib]Note:Algorithms – 矩阵乘法

1.加速递推

2.计算几何

3.……以及其他

[Dr.Lib]Note:Algorithms – LLRB——红黑树的现代实现

一、本文内容

二、算法应用

三、2-3-4树

四、红黑树

1、定义

2、插入算法

3、删除算法

[Dr.Lib]Note:Algorithms – 四边形不等式

2D/1D:d[i,j]=opt{d[i,k-1]+d[k,j]}+w[i,j] (1 \(\leq\) i < k \(\leq\) j \(\leq\) n)

1D/1D:f[i]=opt{f[j]+w[j,i]} (b[i] \(\leq\) j \(\leq\) i-1)

朋友们

TAG CLOUD