学习记录-有后效性的DP状态转移方程（以CH5103和Codeforces24D为例）

前言：

构成动态规划的三要素是 “状态”，“阶段”，“决策”。而动态规划之所以能够从前往后递推，是因为动态规划符合的三个性质：“子问题重叠性”，“无后效性”，“最优子结构性质”。
如果我们已经有了状态的表示和状态转移方程，但是发现dp不满足“无后效性”这一个性质——即部分状态相互联系，相互转移形成了环形，无法确定一个合适的dp阶段，从而沿着某个方向递推。这时，不能再继续按照原来的递推式递推了。

有后效性的解题思路：

高斯消元：https://blog.csdn.net/qq_42785590/article/details/95228501
可以把动态规划的各个状态看做是未知量，状态的转移看做是若干个方程，如果仅仅是“无后效性”这一点不能得到满足，并且状态转移方程都是一次方程，那么我们就可以不进行线性递推，而是用高斯消元直接求出状态转移方程的解。
在题目中，动态规划的转移阶段如果有后效性，那么就是“分阶段带环” 的，我们可以用高斯消元和动态规划相结合，在总体方面采用动态规划，局部带环部分用高斯消元解出互相影响的部分。
下面分别给出两个例子，一个带没有后效性，一个具有后效性。

一个没有后效性的例子：

例题：CH5103
题意：一个方格，每一个格子具有一个权值。从走上角走到右下角，每次只能向下走或者向右走，并且不能超出界限。选择两条路线，求经过的权值和最大的方案。这两条路线都经过的方格只能算一次权值。
n,m<=50。
思路
每一次只能向右走或者向下走，所以可以排除后效性。用d[i][j][k]表示目前两条路线都走了i步，其中第一条路线位于第j行，第二条路线位于第j行的时候，最大权值。那么就有两种情况的递推。
1：如果目前的（j!=k）那么就说明两条路线没有走一个格子上面，就有;
d[i][j][k]=max(max(d[i-1][j-1][k-1],d[i][j][k]),max(d[i-1][j][k-1],d[i-1][j-1][k]))+a[j][i+2-j]+a[i+2-k][k];
2：如果目前的（j==k）那么说明两条路线位于同一个个字上面，那么就有：
d[i][j][k]=max(max(d[i-1][j-1][k-1],d[i][j][k]),max(d[i-1][j][k-1],d[i-1][j-1][k]))+a[j][i+2-j];

代码如下

//#include<bits/stdc++.h>
#include<algorithm>
#include<complex>
#include<iostream>
#include<iomanip>
#include<ostream>
#include<cstring>
#include<string.h>
#include<string>
#include<cstdio>
#include<cctype>
#include<vector>
#include<cmath>
#include<queue>
#include<set>
#include<stack>
#include<map>
#include<cstdlib>
#include<time.h>
#include<ctime>
#include<bitset>
#define pb push_back
#define _fileout freopen("out.txt","w",stdout)
#define _filein freopen("in.txt","r",stdin)
#define ok(i) printf("ok%d\n",i)
using namespace std;
typedef double db;
typedef long long ll;
const double PI = acos(-1.0);
const ll MOD=1e9+7;
const ll NEG=1e9+6;
const int MAXN=2e3+10;
const int INF=0x3f3f3f3f;
const ll ll_INF=9223372036854775807;
const double eps=1e-9;
ll qm(ll a,ll b){ll ret = 1;while(b){if (b&1)ret=ret*a%MOD;a=a*a%MOD;b>>=1;}return ret;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return (a*b)/gcd(a,b);}
int n,m;
int a[50+5][50+5];
int d[50*2+5][50+5][50+5];//表示目前两个路径都走了i步，第一个路径的行数为j，另外一个路径的行数为k时，二者走的最长路径。
int main()
{
    scanf("%d%d",&n,&m);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=m;j++)
        {
            scanf("%d",&a[i][j]);
        }
    }
    d[0][1][1]=a[1][1];
    for(int i=1;i<=m+n-2;i++)
    {
        for(int j=1;j<=min(i+1,n);j++)
        {
            for(int k=1;k<=min(i+1,n);k++)
            {
                int y1=i+2-j;
                int y2=i+2-k;
                if(j==k)d[i][j][k]=max(max(d[i-1][j-1][k-1],d[i-1][j-1][k]),max(d[i-1][j][k-1],d[i-1][j][k]))+a[j][y1];
                else d[i][j][k]=max(max(d[i-1][j-1][k-1],d[i-1][j-1][k]),max(d[i-1][j][k-1],d[i-1][j][k]))+a[j][y1]+a[k][y2];
            }
        }
    }
    printf("%d",d[m+n-2][n][n]);
    return 0;
}

有后效性的例子：

codeforces Problem 24D
题目链接：https://codeforces.com/problemset/problem/24/D
题意
给一个n*m的矩阵，给定起始点x和y，一个人从这个位置出发，每一步可以等概率的向左、向右、向下或者原地不动。求这个人第一次走到最后一行的步数的期望。
n,m<=1000
思路
这道题很显然是一道与概率有关的dp，所以一般来说都是从后往前推。我们用d[i][j]表示目前位于第 i 行第 j 列这个位置，这个人走到最后一行的步数的数学期望。那么可以写出下列递推公式：
j= =1: d[i][j] = 1/3( d[i][j] + d[i][j+1] + d[i+1][j ])+1 ；
j= =m: d[i][j] = 1/3 ( d[i][j] + d[i][j-1] + d[i+1][j] )+1;
1<j<m: d[i][j] = 1/4 ( d[i][j] + d[i][j-1] + d[i][j+1] + d[i+1][j] )+1;
很明显相邻两行之间是没有后效性的，但是同一行之间的状态是相互影响的，所以可以考虑，依旧是从最后一行开始往上一行递推，但是同一行之间可以用多项式来求。利用高斯消元。
其中我们每一次把d[i+1][j]看做是一个常数（因为当我们推到目前这一行的时候已经知道下一行的数据了），根据同一行中每一个d[i][j]的递推公式，我们可以列出一个n维的系数矩阵。利用高斯消元进行多项式求解。
但是高斯消元是O（N ^2）的算法。那么总的算法O（N ^3）的就会T了。但是我们可以观察到系数矩阵的第一行和最后一行只有两个数不是0，即-2/3 和 1/3，第二行到第(m-1)行只有三个数不为0，依次为 1/4, -3/4, 1/4。所以我们高斯消元每行只需要对两个或者三次操作，所以高斯消元的复杂度是O(N)的。所以总的时间复杂度就是(O ^2)。
AC代码

//#include<bits/stdc++.h>
#include<algorithm>
#include<complex>
#include<iostream>
#include<iomanip>
#include<ostream>
#include<cstring>
#include<string.h>
#include<string>
#include<cstdio>
#include<cctype>
#include<vector>
#include<cmath>
#include<queue>
#include<set>
#include<stack>
#include<map>
#include<cstdlib>
#include<time.h>
#include<ctime>
#include<bitset>
#define pb push_back
#define _fileout freopen("out.txt","w",stdout)
#define _filein freopen("in.txt","r",stdin)
#define ok(i) printf("ok%d\n",i)
using namespace std;
typedef double db;
typedef long long ll;
const double PI = acos(-1.0);
const ll MOD=1e9+7;
const ll NEG=1e9+6;
const int MAXN=1e3+10;
const int INF=0x3f3f3f3f;
const ll ll_INF=9223372036854775807;
const double eps=1e-9;
ll qm(ll a,ll b){ll ret = 1;while(b){if (b&1)ret=ret*a%MOD;a=a*a%MOD;b>>=1;}return ret;}
ll gcd(ll a,ll b){return b?gcd(b,a%b):a;}
ll lcm(ll a,ll b){return (a*b)/gcd(a,b);}
int n,m;
int xx,yy;
double d[MAXN][MAXN];
double a[MAXN][MAXN];
double x[MAXN];
int Gauss()
{
    for (int i=1;i<m;i++)//消去未知数
    {
        double rate=a[i+1][i]/a[i][i];
        a[i+1][i]-=rate*a[i][i];
        a[i+1][i+1]-=rate*a[i][i+1];
        x[i+1]-=rate*x[i];
    }   
    x[m]=x[m]/a[m][m]; //最后一行未知数只剩下一个，直接可以得出最后一个未知量的值
    for (int i=m-1;i>0;i--)
    {
        x[i]=(x[i]-x[i+1]*a[i][i+1])/a[i][i];//倒回来求每一行的值
    }
    return 1;
}
int main()
{

    scanf("%d%d%d%d",&n,&m,&xx,&yy);
    if(m==1)
    {
        for(int i=n-1;i>=xx;i--)
        {
            d[i][m]=d[i+1][m]+2;
        }
        printf("%.10f\n",d[xx][yy]);
        return 0;
    }
    for(int i=1;i<=n;i++)
        for(int j=1;j<=m;j++)
            d[i][j]=0.0;
    for(int i=n-1;i>=xx;i--)
    {
        // memcpy(a,aa,sizeof a);
        a[1][1]=(double)-2/3;
        a[1][2]=(double)1/3;
        a[m][m-1]=(double)1/3;
        a[m][m]=(double)-2/3;
        for(int i=2;i<m;i++)
        {
            a[i][i-1]=(double)1/4;
            a[i][i]=(double)-3/4;

            // printf("%.2f\n",a[i][i]);
            a[i][i+1]=(double)1/4;
        }
        x[1]=((double)-1/3)*d[i+1][1]-1;
        x[m]=((double)-1/3)*d[i+1][m]-1;
        for(int j=2;j<m;j++)
        {
            x[j]=((double)-1/4)*d[i+1][j]-1;
        }
        if(Gauss())
        for(int j=1;j<=m;j++)
        {
            d[i][j]=x[j];
            // printf("d[%d][%d]=%.2f ",i,j,d[i][j]);
        }
        // printf("\n");
    }
    printf("%.10f\n",d[xx][yy]);
    return 0;
}