分数规划及斜率优化

$$ A(vec{j})+B(vec{j})F(vec{i})=G(vec{i}) $$

$A$ 和 $B$ 是关于 $vec{j}$ 的函数.

注意,如果确定了 $vec{j}$ ,那么 $G(vec{i})$ 便是一条在平面 $(F(vec{i}),G(vec{i}))$ 上的直线,即一个以 $vec{i}$ 为参数的参数方程

$$ egin{cases}
& y=G(vec{i})=A_j+B_j F(vec{i}) \
& x=F(vec{i})
end{cases} $$

分数规划

考虑

$$ Minimize quad frac{A(vec{x})}{B(vec{x})}=lambda ,; vec{x} in S$$

其中 $vec{x}$ 叫做"解向量",可以理解为一个参数表 $(x_1,x_2,...,x_n)$ ,S是"可行解集合".

$A(vec{x})$ 和 $B(vec{x})$ 是关于 $vec{x}$ 的连续实值函数.

然后我们定义一个函数 $$ P(lambda)=min_{vec{x}in S} left({-B(vec{x})lambda+A(vec{x})} ight) $$

考虑当给定了一个 $lambda$ 后,如何暴力求 $P(lambda)$ .

很明显,直接在集合 $S$ 中枚举 $vec{x}$ ,算出 $B(vec{x})$ 与 $A(vec{x})$ ,算出 $-B(vec{x})lambda+A(vec{x})$ ,取最小值即可.

然后我们注意看这个过程. 我们枚举了直线 $$y=-B(vec{x})lambda+A(vec{x})$$

其中自变量 $lambda$ ,因变量 $y$ .并计算出它的 $y$ 值.

这相当于我们使用一条垂直于水平坐标轴($lambda$轴)的直线,与所有直线 $y=-B(vec{x})lambda+A(vec{x})$ 取交点,然后取得交点的纵坐标最小的那个.

回到原问题,求 $lambda$ 的最小值.

我们发现,对于一条直线 $y=-B(vec{x})lambda+A(vec{x})$ ,它的解(使用可行解 $vec{x}$ 所能给出的 $lambda$ )正好对应了直线在 $lambda$ 轴上的截距.

因此,原问题等价于:

$$找出一个vec{x},使得(由vec{x}确定的)直线y=-B(vec{x})lambda+A(vec{x})在lambda轴上的截距最小.$$

或者有时原问题要求的是最小值而不求具体方案,则原问题等价于

$$找出一个lambda_0 ,使得存在一个可行解vec{x}in S,\ 直线 y=-B(vec{x})lambda+A(vec{x}) 在lambda轴上的截距为lambda_0,\ 且不存在其它可行解vec{x}in S,使得直线截距比lambda_0 小.$$

这个时候就要用到一些既定的性质了.

一般分数规划题目会有$$forall _{vec{x}in S} ; B(vec{x})>0$$

于是,我们发现直线 $y=-B(vec{x})lambda+A(vec{x})$ 的斜率总是比0小的.

于是 $P(lambda)$ 是单调的.

那么我们的目标就很明确了:

$$求出函数P(lambda)与lambda 轴的左交点.$$

红线就是 $lambda$ 函数图像.可以看到所有直线的$lambda$轴交点都在右边.

单调就可以二分查找$lambda$啦.......

但是如何计算$P(lambda)$呢? 刚才的枚举完全不可做啊.....

呃,这就要具体问题具体分析了.......

一般都会和原算法扯上关系,比如最优比率路径使用最短路径算法计算$P(lambda)$,最优比率生成树使用最小生成树算法计算$P(lambda)$,最优比率割使用最小割计算$P(lambda)$.......

AC POJ 2728 最优比率生成树分数规划+Prim

  1 #include <cstdio>
  2 #include <fstream>
  3 #include <iostream>
  4  
  5 #include <cstdlib>
  6 #include <cstring>
  7 #include <algorithm>
  8 #include <cmath>
  9  
 10 #include <queue>
 11 #include <vector>
 12 #include <map>
 13 #include <set>
 14 #include <stack>
 15 #include <list>
 16  
 17 typedef unsigned int uint;
 18 typedef long long int ll;
 19 typedef unsigned long long int ull;
 20 typedef double db;
 21  
 22 using namespace std;
 23  
 24 inline int getint()
 25 {
 26     int res=0;
 27     char c=getchar();
 28     bool mi=false;
 29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 31     return mi ? -res : res;
 32 }
 33 inline ll getll()
 34 {
 35     ll res=0;
 36     char c=getchar();
 37     bool mi=false;
 38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 40     return mi ? -res : res;
 41 }
 42  
 43 db eps=1e-80;
 44 inline bool feq(db a,db b)
 45 { return fabs(a-b)<eps; }
 46 
 47 template<typename Type>
 48 inline Type avg(const Type a,const Type b)
 49 { return a+((b-a)/2); }
 50 
 51 //===================================================================
 52 //===================================================================
 53 //===================================================================
 54 //===================================================================
 55 
 56 const int INF=(1<<30)-1;
 57 
 58 int n;
 59 
 60 db M[1050][1050]; //adjust table.
 61 db dist[1050][1050];
 62 
 63 int px[1050];
 64 int py[1050];
 65 int ph[1050];
 66 
 67 db Dist(db ax,db ay,db bx,db by)
 68 { return sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by)); }
 69 
 70 bool intree[1050];
 71 db dst[1050];
 72 
 73 db MST(db f)
 74 {
 75     for(int i=0;i<n;i++)
 76     for(int j=i;j<n;j++)
 77     M[i][j]=M[j][i]=fabs(ph[i]-ph[j])-f*dist[i][j];
 78     
 79     for(int i=0;i<n;i++)
 80     dst[i]=1e40,intree[i]=false;
 81     
 82     db res=0.0;
 83     
 84     intree[0]=true;
 85     for(int i=0;i<n;i++)
 86     dst[i]=min(dst[i],M[0][i]);
 87     
 88     for(int d=1;d<n;d++)
 89     {
 90         db mi=1e40;
 91         int mp=-1;
 92         for(int i=0;i<n;i++)
 93         if(!intree[i] && mi>dst[i])
 94         {
 95             mi=dst[i];
 96             mp=i;
 97         }
 98         
 99         intree[mp]=true;
100         res+=dst[mp];
101         
102         for(int i=0;i<n;i++)
103         if(!intree[i])
104         dst[i]=min(dst[i],M[mp][i]);
105     }
106     
107     return res;
108 }
109 
110 int main()
111 {
112     eps=1e-4;
113     
114     while(true)
115     {
116         n=getint();
117         if(n==0) break;
118         
119         for(int i=0;i<n;i++)
120         for(int j=0;j<n;j++)
121         M[i][j]=INF;
122         
123         for(int i=0;i<n;i++)
124         {
125             px[i]=getint()-1;
126             py[i]=getint()-1;
127             ph[i]=getint();
128         }
129         
130         for(int i=0;i<n;i++)
131         for(int j=i;j<n;j++)
132         dist[i][j]=dist[j][i]=Dist(px[i],py[i],px[j],py[j]);
133         
134         db l=0,r=1e7;
135         while(fabs(r-l)>eps)
136         {
137             db mid=(l+r)*0.5;
138             if(MST(mid)>0.0) l=mid;
139             else r=mid;
140         }
141         
142         printf("%.3f
",l);
143     }
144     
145     return 0;
146 }

View Code

具体用分数规划解决问题的时候就直接

1.确定比式 $frac{A(vec{x})}{B(vec{x})}=lambda$

2.构造函数 $P(lambda)=min(A(vec{x})+B(vec{x})lambda )$.

这道题是 $P(lambda)=min(vec{a}cdotvec{x}-vec{b}cdotvec{x} lambda)$ ,

其中 $vec{a}$ 代表高度差数组(答案的分子部分), $vec{b}$ 代表距离数组(答案的分母部分).

3.整理表达式.这道题是 $min((vec{a}-lambda vec{b})cdotvec{x})$

4.考虑用某种常规方法拿到 $P(lambda)$. 这道题由 $min((vec{a}-lambda vec{b})cdotvec{x})$ 的意义,知这个值可以用MST算法求出.

5.二分/三分 $lambda$ 达到精度要求后直接输出.

斜率优化

一般用在DP中,求$$min_{0le j < i} A(j)+B(j)cdot F(i) $$

也可以求 $max$ ,这里只讨论 $min$ .

注意到每一个j确定了一条直线 $ y=B(j)x+A(j) $

当我们用暴力枚举 $j$ 计算状态 $i$ 的转移的时候,实际上我们在枚举直线 $y=B(j)x+A(j)$ ,然后把 $x=F(i)$ 代入,求得y的最小值..

根据上边某个图可以知道其实就是求一个下凸壳(或上凸壳)与直线 $x=F(i)$ 的交点.

这样我们维护凸壳就行了$=omega=$(说得倒简单!)

呃嗯,斜率优化是这样一种东西..

我们需要找一些性质......

我们就可以比较方便地解决问题.....

1.斜率有单调性.

此时所有的直线的斜率 $B(j)$ 按照 $j$ 的顺序单调递减(或递增,这里只讨论递减). 画一下图就知道,我们需要维护上凸壳(下半平面交),并且由于斜率递减的原因,我们总是从后往前(从右往左)剔除一些没用的直线.

具体的,考虑新加入一条直线 $l_0$ ,而最右边的直线,按照斜率记为 $l_1$ , $l_2$ ,($l_1$ 斜率小于 $l_2$ 斜率,即 $l_1$ 是半平面交的直线栈中最后一个元素 ),然后判断 $l_0$ 与 $l_2$ 的交点是不是在 $l_1$ 的正下方.

如果是,那么 $l_1$ 一定不会被选中成为最优的转移状态,因为 $l_1$ 与 $l_2$ 已经完全把它的正下方遮住了.

重复这个剔除直线的操作直到不能剔除( $l_0$ 与 $l_2$ 交点在 $l_1$ 上方 )为止.

这样我们就可以直接像GRAHAM一样维护了......

配合二分(三分)查找(就是二分(三分)找到一条直线 $l$ ,处在它左边的(即凸壳栈中和它左相邻的)直线 $l^{-}$ 和 $x=F(i)$ 的交点在 $l$ 上方,处在它右边的直线 $l^{+}$ 和 $x=F(i)$ 的交点也在 $l$ 上方. ),可以得到它的值. 复杂度 $O(nlogn)$

2.决策参数 $F(i)$ 有单调性

这个东西需要斜率单调.

斜率单调以后,如果 $F(i)$ 随着 $i$ 表现单调的话,我们可以从凸壳某一侧排除一些不再能够成为解的直线.

如果凸壳最左边的直线 $l_0$ 比右边的直线 $l_1$ 的解更差,那么我们可以断定,对于以后的 $x=F(i)$ ,这条直线都不会成为解. 于是可以搞一个单调队列,相当于指针扫描,队列在不断插入新边的同时也不断地剔除左边没用的直线.可以在 $O(n)$ 时间解决问题.

如果两个条件都不满足,也能做,因为只要是类似于 $min$ 或 $max$ 函数中的表达式是直线形式,都可以形成凸壳,这样我们需要splay/线段树维护凸壳,二分查找.

有些神题还要可持久化.......

AC BZOJ 1010

  1 #include <cstdio>
  2 #include <fstream>
  3 #include <iostream>
  4  
  5 #include <cstdlib>
  6 #include <cstring>
  7 #include <algorithm>
  8 #include <cmath>
  9  
 10 #include <queue>
 11 #include <vector>
 12 #include <map>
 13 #include <set>
 14 #include <stack>
 15 #include <list>
 16  
 17 typedef unsigned int uint;
 18 typedef long long int ll;
 19 typedef unsigned long long int ull;
 20 typedef double db;
 21  
 22 using namespace std;
 23  
 24 inline int getint()
 25 {
 26     int res=0;
 27     char c=getchar();
 28     bool mi=false;
 29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 31     return mi ? -res : res;
 32 }
 33 inline ll getll()
 34 {
 35     ll res=0;
 36     char c=getchar();
 37     bool mi=false;
 38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 40     return mi ? -res : res;
 41 }
 42  
 43 db eps=1e-20;
 44 inline bool feq(db a,db b)
 45 { return fabs(a-b)<eps; }
 46 
 47 template<typename Type>
 48 inline Type avg(const Type a,const Type b)
 49 { return a+((b-a)/2); }
 50 
 51 //==========================================================
 52 //==========================================================
 53 //==========================================================
 54 //==========================================================
 55 
 56 ll INF=((ll)1<<56)-1;
 57 
 58 ll n,L;
 59 
 60 ll d[50050];
 61 ll c[50050];
 62 ll s[50050];
 63 
 64 ll bd[50050];
 65 
 66 ll F[50050];
 67 ll A[50050];
 68 ll B[50050];
 69 ll K[50050];
 70 
 71 ll C;
 72 
 73 /*
 74 inline ll F(int i)
 75 { return i+s[i]; }
 76 
 77 inline ll A(int i)
 78 { return F(i)*F(i)-2*C*F(i); }
 79 
 80 inline ll B(int i)
 81 { return F(i)*F(i)+2*C*F(i)+d[i]; }
 82 
 83 inline ll K(int i)
 84 { return -2*F(i); }
 85 */
 86 
 87 int stk[50050],hp=0,tp=0;
 88 inline int getans(int i)
 89 {
 90     while( hp+1<tp && B[stk[hp]]-B[stk[hp+1]] >= (K[stk[hp+1]]-K[stk[hp]])*F[i] )
 91     hp++;
 92     return stk[hp];
 93 }
 94 
 95 inline void addline(ll b,ll k,int i)
 96 {
 97     while(  hp+1<tp &&
 98             (b-B[stk[tp-1]])*(k-K[stk[tp-2]]) >=
 99             (b-B[stk[tp-2]])*(k-K[stk[tp-1]]) )
100         tp--;
101     stk[tp++]=i;
102 }
103 
104 
105 
106 int main()
107 {
108     n=getint();
109     L=getint();
110     C=L+1;
111     
112     for(int i=1;i<=n;i++) c[i]=getint();
113     for(int i=1;i<=n;i++) s[i]=s[i-1]+c[i];
114     
115     tp=1;
116     F[0]=0+s[0];
117     A[0]=(F[0]-(C<<1))*F[0];
118     K[0]=-(F[0]<<1);
119     B[0]=(F[0]+(C<<1))*F[0]+d[0];
120         
121     for(int i=1;i<=n;i++)
122     {
123         F[i]=i+s[i];
124         A[i]=(F[i]-(C<<1))*F[i];
125         K[i]=-(F[i]<<1);
126         
127         int j=getans(i);
128         d[i]=B[j]+K[j]*F[i];
129         d[i]+=A[i]+C*C;
130         
131         B[i]=(F[i]+(C<<1))*F[i]+d[i];
132         addline(B[i],K[i],i);
133     }
134     
135     printf("%lld
",d[n]);
136     
137     return 0;
138 }

View Code

裸的斜率优化.

1.注意队列的端界. 我们允许队列中只存在一个元素(即,凸壳中有且仅有一条直线). 即便当我们加入了另一条直线,使凸壳绝不可能只含一条直线,我们也需要那样判断以免出现坑爹情况,尤其是使用单调队列的时候.

2.公式千万别推错啊!.....

3.注意下标之间的关系......栈中存的是状态编号.........

AC BZOJ 3156

  1 #include <cstdio>
  2 #include <fstream>
  3 #include <iostream>
  4  
  5 #include <cstdlib>
  6 #include <cstring>
  7 #include <algorithm>
  8 #include <cmath>
  9  
 10 #include <queue>
 11 #include <vector>
 12 #include <map>
 13 #include <set>
 14 #include <stack>
 15 #include <list>
 16  
 17 typedef unsigned int uint;
 18 typedef long long int ll;
 19 typedef unsigned long long int ull;
 20 typedef double db;
 21  
 22 using namespace std;
 23  
 24 inline int getint()
 25 {
 26     int res=0;
 27     char c=getchar();
 28     bool mi=false;
 29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 31     return mi ? -res : res;
 32 }
 33 inline ll getll()
 34 {
 35     ll res=0;
 36     char c=getchar();
 37     bool mi=false;
 38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 40     return mi ? -res : res;
 41 }
 42  
 43 db eps=1e-80;
 44 inline bool feq(db a,db b)
 45 { return fabs(a-b)<eps; }
 46 
 47 template<typename Type>
 48 inline Type avg(const Type a,const Type b)
 49 { return a+((b-a)/2); }
 50 
 51 //===================================================================
 52 //===================================================================
 53 //===================================================================
 54 //===================================================================
 55 
 56 int n;
 57 ll a[1005000];
 58 ll d[1005000];
 59 
 60 ll A[1005000];
 61 ll K[1005000];
 62 
 63 int stk[1005000],hp,tp;
 64 int getans(ll v)
 65 {
 66     while( hp+1<tp &&
 67             A[stk[hp]]-A[stk[hp+1]] >= (K[stk[hp+1]]-K[stk[hp]])*v )
 68         hp++;
 69     return stk[hp];
 70 }
 71 void pushline(ll a,ll k,int i)
 72 {
 73     while( hp+1<tp &&
 74             (a-A[stk[tp-1]])*(k-K[stk[tp-2]]) >= (a-A[stk[tp-2]])*(k-K[stk[tp-1]]) )
 75         tp--;
 76     stk[tp++]=i;
 77 }
 78 
 79 
 80 int main()
 81 {
 82     n=getint();
 83     for(int i=1;i<=n;i++) a[i]=getint();
 84     
 85     hp=0; tp=1;
 86     
 87     for(int i=1;i<=n;i++)
 88     {
 89         int j=getans(i);
 90         d[i]=A[j]+K[j]*i;
 91         d[i]+=a[i]+(ll)(i-1)*i/2;
 92         
 93         A[i]=d[i]+(ll)(i+1)*i/2;
 94         K[i]=-i;
 95         
 96         pushline(A[i],K[i],i);
 97     }
 98 
 99     printf("%lld
",d[n]);
100     return 0;
101 }

View Code

由于循环变量用的int,结果乘起来爆精度了.......WA了好几发.......

为什么都这时候了还错在这种sb问题上........

而且推式子的时候居然傻逼地把min提一个1/2常数出来结果把常数项从min里拿到外面又忘记乘.........妈呀.......

AC BZOJ 3675 APIO2014 T2

比较裸的斜率优化.

  1 #include <cstdio>
  2 #include <fstream>
  3 #include <iostream>
  4  
  5 #include <cstdlib>
  6 #include <cstring>
  7 #include <algorithm>
  8 #include <cmath>
  9  
 10 #include <queue>
 11 #include <vector>
 12 #include <map>
 13 #include <set>
 14 #include <stack>
 15 #include <list>
 16  
 17 typedef unsigned int uint;
 18 typedef long long int ll;
 19 typedef unsigned long long int ull;
 20 typedef double db;
 21  
 22 using namespace std;
 23  
 24 inline int getint()
 25 {
 26     int res=0;
 27     char c=getchar();
 28     bool mi=false;
 29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 31     return mi ? -res : res;
 32 }
 33 inline ll getll()
 34 {
 35     ll res=0;
 36     char c=getchar();
 37     bool mi=false;
 38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
 39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
 40     return mi ? -res : res;
 41 }
 42  
 43 db eps=1e-20;
 44 inline bool feq(db a,db b)
 45 { return fabs(a-b)<eps; }
 46 
 47 template<typename Type>
 48 inline Type avg(const Type a,const Type b)
 49 { return a+((b-a)/2); }
 50 
 51 //==============================================================================
 52 //==============================================================================
 53 //==============================================================================
 54 //==============================================================================
 55 
 56 const ll INF=((ll)1<<56)-1;
 57 
 58 int n;
 59 int m;
 60 
 61 int a[105000];
 62 ll s[105000];
 63 ll s2[105000];
 64 ll d[2][105000];
 65 
 66 ll b[2][105000];
 67 
 68 int st[105000];
 69 int h,t;
 70 void INIT(){ h=t=0; }
 71 
 72 inline void Insert(int l,ll B,ll K,int id)
 73 {
 74     while(h<t-1 && (B-b[l][st[t-1]])*(K-s[st[t-2]]) >=
 75                     (B-b[l][st[t-2]])*(K-s[st[t-1]]) ) t--;
 76     st[t++]=id;
 77 }
 78 
 79 inline int Get(int l,ll x)
 80 {
 81     while( h<t-1 && b[l][st[h+1]]-b[l][st[h]] >=
 82                     (s[st[h]]-s[st[h+1]])*x ) h++;
 83     return st[h];
 84 }
 85 
 86 int main()
 87 {
 88     n=getint();
 89     m=getint();
 90     for(int i=1;i<=n;i++)
 91     s[i]=s[i-1]+(a[i]=getint()),s2[i]=s[i]*s[i];
 92     
 93     int cur=0;
 94     
 95     for(int i=0;i<=n;i++)
 96     b[1][i]=-s2[i];
 97     
 98     memset(d[!cur],0,sizeof(ll)*(n+1));
 99     
100     for(int r=1;r<=m;r++)
101     {    
102         memset(d[cur],0,sizeof(ll)*(n+1));
103         INIT();
104         Insert(!cur,b[!cur][0],s[0],0);
105         
106         for(int i=1;i<=n;i++)
107         {
108             int v=Get(!cur,s[i]);
109             d[cur][i]=d[!cur][v]-s2[v]+s[v]*s[i];
110             
111             b[cur][i]=d[cur][i]-s2[i];
112             
113             Insert(!cur,b[!cur][i],s[i],i);
114         }
115         
116         cur=!cur;
117     }
118     printf("%lld
",d[!cur][n]);
119     
120     return 0;
121 }

View Code

唯一的不同在于多了一维....

状态是从d[cur-1]转移至d[cur],所以单调队列里的所有东西都是d[cur-1]里的内容.

然后为了节省空间.....我这里写了滚动数组.......

注意倒式子的时候的符号...... "不断论证推导的每一步是否正确,为什么正确."

注意这道题要用到一个性质:一旦确定了所有切割点,那么切割顺序不影响最终结果. 如何证明?

AC BZOJ1911

 1 #include <cstdio>
 2 #include <fstream>
 3 #include <iostream>
 4  
 5 #include <cstdlib>
 6 #include <cstring>
 7 #include <algorithm>
 8 #include <cmath>
 9  
10 #include <queue>
11 #include <vector>
12 #include <map>
13 #include <set>
14 #include <stack>
15 #include <list>
16  
17 typedef unsigned int uint;
18 typedef long long int ll;
19 typedef unsigned long long int ull;
20 typedef double db;
21  
22 using namespace std;
23  
24 inline int getint()
25 {
26     int res=0;
27     char c=getchar();
28     bool mi=false;
29     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
30     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
31     return mi ? -res : res;
32 }
33 inline ll getll()
34 {
35     ll res=0;
36     char c=getchar();
37     bool mi=false;
38     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
39     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
40     return mi ? -res : res;
41 }
42 
43 //==============================================================================
44 //==============================================================================
45 //==============================================================================
46 //==============================================================================
47 
48 int n;
49 ll A,B,C;
50 
51 ll a[1005000];
52 ll s[1005000];
53 ll d[1005000];
54 
55 ll stb[1005000];
56 ll stk[1005000];
57 int sh,st;
58 void Insert(ll b,ll k)
59 {
60     while(st>1 && (stk[st-2]-k)*(stb[st-2]-stb[st-1])<(stk[st-2]-stk[st-1])*(stb[st-2]-b)) st--;
61     stb[st]=b;
62     stk[st]=k;
63     st++;
64 }
65 
66 ll Get(ll v)
67 {
68     while(sh<st-1 && stb[sh]+stk[sh]*v < stb[sh+1]+stk[sh+1]*v) sh++;
69     return stb[sh]+stk[sh]*v;
70 }
71 
72 int main()
73 {
74     n=getint();
75     A=getint();
76     B=getint();
77     C=getint();
78     
79     for(int i=1;i<=n;i++)
80     a[i]=getint();
81     
82     for(int i=1;i<=n;i++)
83     s[i]=s[i-1]+a[i];
84     
85     Insert(0,0);
86     for(int i=1;i<=n;i++)
87     {
88         d[i]=Get(s[i])+A*s[i]*s[i]+B*s[i]+C;
89         Insert(A*s[i]*s[i]-B*s[i]+d[i],-2*A*s[i]);
90     }
91     
92     printf("%lld
",d[n]);
93     
94     return 0;
95 }

View Code

推式子整个推错了居然没发现QAQ

注意正负号,用直线式算出来那个坐标值一般都带负号的.

急着A题结果把减号打成乘号居然也没有发现啊QAQ

打题的时候一定要冷静....冷静........

AC BZOJ1096 裸的斜率优化

 1 #include <cstdio>
 2 #include <fstream>
 3 #include <iostream>
 4  
 5 #include <cstdlib>
 6 #include <cstring>
 7 #include <algorithm>
 8 #include <cmath>
 9  
10 #include <queue>
11 #include <vector>
12 #include <map>
13 #include <set>
14 #include <stack>
15 #include <list>
16  
17 typedef unsigned int uint;
18 typedef long long int ll;
19 typedef unsigned long long int ull;
20 typedef double db;
21 typedef long double ldb;
22  
23 using namespace std;
24  
25 inline int getint()
26 {
27     int res=0;
28     char c=getchar();
29     bool mi=false;
30     while((c<'0' || c>'9')/* && !feof(stdin)*/) mi=(c=='-'),c=getchar();
31     while('0'<=c && c<='9'/* && !feof(stdin)*/) res=res*10+c-'0',c=getchar();
32     return mi ? -res : res;
33 }
34 inline ll getll()
35 {
36     ll res=0;
37     char c=getchar();
38     bool mi=false;
39     while(c<'0' || c>'9') mi=(c=='-'),c=getchar();
40     while('0'<=c && c<='9') res=res*10+c-'0',c=getchar();
41     return mi ? -res : res;
42 }
43 
44 //==============================================================================
45 //==============================================================================
46 //==============================================================================
47 //==============================================================================
48 
49 
50 
51 struct factory
52 { ll x,p,c; }F[1005000];
53 bool fcmp(const factory&a,const factory&b)
54 { return a.x<b.x; }
55 ll s[1005000];
56 ll t[1005000];
57 int n;
58 
59 ll B[1005000];
60 ll K[1005000];
61 int qh,qt;
62 ll Get(ll x)
63 {
64     while(qh<qt-1 && B[qh]-B[qh+1] > x*(K[qh+1]-K[qh]) ) qh++;
65     return B[qh]+K[qh]*x;
66 }
67 ll Insert(ll b,ll k)
68 {
69     while( qt>qh+1 && 
70           (K[qt-2]-k)*(B[qt-2]-B[qt-1])<=(K[qt-2]-K[qt-1])*(B[qt-2]-b) ) qt--;
71     B[qt]=b; K[qt]=k; qt++; 
72 }
73 
74 ll d[1005000];
75 
76 int main()
77 {
78     n=getint();
79     for(int i=1;i<=n;i++)
80     F[i].x=getint(),F[i].p=getint(),F[i].c=getint();
81     sort(F,F+n,fcmp);
82     for(int i=1;i<=n;i++) s[i]=s[i-1]+F[i].p;
83     for(int i=1;i<=n;i++) t[i]=t[i-1]+F[i].x*F[i].p;
84     
85     Insert(0,0);
86     for(int i=1;i<=n;i++)
87     {
88         ll v=Get(F[i].x);
89         d[i]=v+F[i].x*s[i]-t[i]+F[i].c;
90         Insert(d[i]+t[i],-s[i]);
91     }
92     
93     printf("%lld
",d[n]);
94     
95     return 0;
96 }

View Code

又被 long long 坑掉 TAT....

随手在结构体里打了int结果就WA了一发.....

另一种斜率优化

考虑方程:$$ max_{0le j<i}{ x_j A_i + Y_j B_i }$$

此时我们多了一个与i相关的控制变量,所以用于评判结果的直线不再是垂直于x轴的线.

考虑平面上的点 $(X_j,Y_j)$ . 对于一个确定的 $A_i$ 以及 $B_i$ , 我们可以构造一条直线 $ p=x A_i + y B_i $ .

这样,如果直线经过了一个点 $(X_j,Y_j)$ ,这条直线的 $p$ 值就等于我们所要求的 $ X_j A_i + Y_j B_i $ .

于是,我们可以想象一下,一条直线 $p=A_i x + N_i y$ ,一开始p为一个非常大的数,然后p一直减少,这样直线就会一直向下平移,直到我们碰到一个点 $(X_j,Y_j)$ ,那么我们就取得了最大的那个 $p$ ,就是我们要求的结果.

但是如何计算呢? 我们可以任选一个p的值,搞一条直线出来. 考虑点到直线的距离

$$ D = frac{| A_i X_j + B_i Y_j + p |}{sqrt{A_i^2 + B_i^2} } $$

如果把绝对值去掉,就可以根据符号判断点在直线的上方还是下方.

进一步,由于 $sqrt{A_i^2+B_i^2}$ 不影响点到直线距离的相对大小和符号(我们把所有点到直线的距离同时乘上这个数),得到一个可行判据

$$ D^prime = A_i X_j + B_i Y_j + p $$

画个图就知道,如果我们按照凸壳上的点从左到右的顺序算这个 $D^prime$ 值,会发现它先增大,后减小,是一个上凸函数.

用三分法找到D值最大的那个点就可以直接算出答案了.

值得注意的是,本分析中使用的是Max函数,维护上凸壳. 而上边另一个方程使用Min函数,维护的也是上凸壳.

两个模型是否能相互转化呢?

原来cnblogs的LaTeX不能使用CopyPaste...

从大仙那里讨教来了一些神奇的东西......

china-pub

同方数据库

柳高也买有数据库....

还有《组合数学》 by richard a.brualdi