bzoj2817[ZJOI2012]波浪 DP+高精度

题意：给你三个数n,m,q。问你一个1-n的排列中，两两差值之和大于等于m的概率是多少，保留q位小数,$q<=30$。
可以发现，30位肯定要高精度。。
然后，我们先把所有的可能计算出来，然后除以n!就可以了。
现在问题是如何计算所有情况的总和，那么我们可以直接dp。
可以发现，如果把排列变成连续一段，可以直接O（1）统计答案。
那么我们设f[i][j][k][s]表示做到第i位，当前贡献为j，有k个连续段，边界是s，因为边界的计算不同于段内，所以需要多余记录。
那么有三种情况：
1.新加入的数字新开一个段，贡献为-2i
2.新加入的数延伸原来的段,贡献为2i
3.新加入的数字把两个段链接在一起，贡献为0
至于为什么请自行思考（随便玩玩就能出来）
那么就只剩下一些细节了。。
30位时直接高精度显然常数爆炸，所以我们可以用float128.
但是全部都用显然会T，那么我们小一点的用double，大一点的用float

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<iostream>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;
const int N=1e5+5;
typedef long long ll;
const int mo=4500,MX=9000;
namespace normal
{
    typedef double db;
    db f[2][101][9005][3];
 } 
namespace extend
{
    typedef __float128 db;
    db f[2][101][9005][3];
}
int sta[N],top;
template <class node>
inline void print(node a,int p,char c=0)
{
    ll s=(ll)a;
    a-=s;
    for(top=0;s;s/=10)sta[++top]=s%10;
    if(!top)sta[++top]=0;
    int x=top;
    fo(i,1,p+1)a*=10,sta[++top]=(int)a,a-=sta[top];
    int g=0;
    if (sta[top]>=5)g=1;
    top--;
    fd(i,top,1)
    {
        sta[i]+=g;
        g=0;
        if (sta[i]>=10)sta[i]-=10,g=1;
     } 
     if (g)
     {
        ++top;++x;
        fd(i,top,2)sta[i]=sta[i-1];
        sta[1]=1;
     }
     fo(i,1,x)putchar(sta[i]+'0');
     if (!p)
     {
        putchar(c);
        return;
     }
     else putchar('.');
     fo(i,x+1,top)putchar(sta[i]+'0');
     putchar(c);
}
int n,m,q;
template <class node>
inline void solve(node f[][101][9005][3])
{
    int x=0,last=0;
    f[x][1][-2+mo][0]=1;
    f[x][1][-1+mo][1]=2;
    f[x][1][0+mo][2]=1;
    fo(i,2,n)
    {
        x=last^1;
        memset(f[x],0,sizeof(f[x]));
        fo(j,0,mo*2)
        {
            fo(k,1,n-1)
            {
                if (f[last][k][j][0])
                {
                    node tmp=f[last][k][j][0];
                    if (j-2*i>=0)f[x][k+1][j-2*i][0]+=tmp*(k+1);
                    if (j+2*i<=MX)f[x][k-1][j+2*i][0]+=tmp*(k-1);
                    f[x][k][j][0]+=tmp*2*k;
                    if (j+i<=MX)f[x][k][j+i][1]+=tmp*2;
                    if (j-i>=0)f[x][k+1][j-i][1]+=tmp*2;
                }
                if (f[last][k][j][1])
                {
                    node tmp=f[last][k][j][1];
                    if (j-2*i>=0)f[x][k+1][j-2*i][1]+=tmp*k;
                    if (j+2*i<=MX)f[x][k-1][j+2*i][1]+=tmp*(k-1);
                    f[x][k][j][1]+=tmp*(2*k-1);
                    if (j+i<=MX)f[x][k][j+i][2]+=tmp;
                    if (j-i>=0)f[x][k+1][j-i][2]+=tmp;
                }
                if (f[last][k][j][2])
                {
                    node tmp=f[last][k][j][2];
                    if (j-2*i>=0)f[x][k+1][j-2*i][2]+=tmp*(k-1);
                    if (j+2*i<=MX)f[x][k-1][j+2*i][2]+=tmp*(k-1);
                    f[x][k][j][2]+=tmp*2*(k-1);
                }
            }
        }
        last=x;
    }
    node res=0;
    fo(i,m+mo,mo*2)res+=f[x][1][i][2];
    fo(i,1,n)res/=(node)i;
    print(res,q,'\n');
}
int main()
{
    freopen("river.in","r",stdin);
    freopen("river.out","w",stdout);
    scanf("%d%d%d",&n,&m,&q);
    if (q<=8)solve(normal::f);
    else solve(extend::f);
    return 0;
}