Codeforces989E——A Trance of Nightfall（矩阵乘法+倍增优化）

传送门

大意： 一群人逛景点，总共有n 个景点，坐标分别为$(x_i,y_i)$，他们每次移动按照下面的顺序操作：
1、选择一条直线，要求直线经过现在的位置和至少两个景点（如果现在在某个景点那 里，也算一个）如果有多条直线满足要求，等概率选择一条。
2、在选择的这条直线中，等概率选择一个直线覆盖了的景点移动过去，如果目前在景 点上，也有可能停住不动。
总共有q次询问，第 i 次询问从一个你选的任意点出发（可以不是景点），然后 连续移动 $m_i$步，最后到达 $t_i$的最大概率是多少。

自己完全不会啊，，，靠着膜zxy神仙才会的。。。。

主要的思路看神仙博客就可以了

我主要处理一些

$f[i][j]$表示走了i步到j这个位置的概率，初始化$f[0][m]=1$
如果从x 有cnt条直线，y所在的直线有num 个点，那么从x 到y 的

概率是$1/cnt/num$

然后我们发现可以通过矩阵乘法优化

比如$i$走一步到$k$的概率乘上$k$走一步到$j$的概率就是从$i$走两步到$j$的概率（当然只是针对这三个点之间方案，当然也可能有其他可能，但没有影响，矩阵乘法都会计算在内的）

但是发现直接乘一次是$n^3logm$的，复杂度会爆炸

所以考虑倍增预处理优化

因为我们不可能直接预处理出1e4以内的所以可能，预处理就要爆炸

所以我们考虑倍增优化

由于矩阵自乘满足交换律

所以我们可以处理出$logm$个矩阵表示$2^i$次方的结果

这样可以在每次logm的结果内找到答案

然后注意有两种情况

如果我们直接从一个景点开始走的话是直接算的

但是如果是从一个不是景点的点开始的话就首先要走一步走到景点，所以还要分类讨论一下

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define eps 1e-9
inline int read(){
    char ch=getchar();
    int res=0,f=1;
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch))res=(res<<3)+(res<<1)+(ch^48),ch=getchar();
    return res*f;
}
const int N=205;
int n,x[N],y[N],cnt[N];
double f[16][N][N],g[N],tmp[N];
bool vis[N];
vector<int> G[N][N];
vector<pair<int,int> >line;
inline bool check(int u,int v,int i){
    return (x[i]-x[v])*(y[v]-y[u])==(x[v]-x[u])*(y[i]-y[v]);
}
int main(){
    n=read();
    for(int i=1;i<=n;++i){
        x[i]=read(),y[i]=read();
    }
    for(int i=1;i<=n;++i){
        memset(vis,0,sizeof(vis));
        for(int j=1;j<=n;j++){
            if(i==j||vis[j]) continue;
            cnt[i]++;
            for(int k=1;k<=n;k++){
                if(check(i,j,k))G[i][j].push_back(k),vis[k]=true;
            }
            line.push_back(make_pair(G[i][j][0],G[i][j][1]));
        }
    }
    sort(line.begin(),line.end());
    line.erase(unique(line.begin(),line.end()),line.end());
    for(int i=0;i<line.size();i++){
        vector<int> vec=G[line[i].first][line[i].second];
        for(int j=0;j<vec.size();j++){
            for(int k=0;k<vec.size();k++){
                f[0][vec[j]][vec[k]]+=1.0/(1.0*vec.size());
            }
        }
    }
    for(int i=1;i<=n;i++){
        for(int j=1;j<=n;j++){
            f[0][i][j]/=cnt[i];
        //    cout<<f[0][i][j]<<" ";
        }
//        puts("");
    }
//    for(int i=1;i<=n;i++)cout<<cnt[i]<<" ";
    for(int i=1;i<=15;i++){
        for(int j=1;j<=n;j++){
            for(int k=1;k<=n;k++){
                if(f[i-1][j][k]>1e-6){
                    for(int p=1;p<=n;p++){
                        f[i][j][p]+=f[i-1][j][k]*f[i-1][k][p];
                    }
                }
            }
        }
    }
    int q=read();
    for(int cas=1;cas<=q;cas++){
        int des=read(),step=read()-1;
        memset(g,0,sizeof(g));
        g[des]=1;
        for(int i=0;i<=15;i++){
            if((1<<i)>step)break;
            if((1<<i)&step){
                memset(tmp,0,sizeof(tmp));
                for(int j=1;j<=n;j++){
                    if(g[j]>eps){
                        for(int k=1;k<=n;k++){
                            tmp[k]+=f[i][k][j]*g[j];
                        }
                    }
                }
                memcpy(g,tmp,sizeof(tmp));
            }
        }
        double ans=0;
        for(int i=0;i<line.size();i++){
            vector<int> vec=G[line[i].first][line[i].second];
            double sum=0;
            for(int j=0;j<vec.size();j++){
                sum+=g[vec[j]];
            }
            sum/=vec.size();
            ans=max(ans,sum);
        }
        memset(tmp,0,sizeof(tmp));
        for(int i=1;i<=n;i++){
            if(g[i]>eps){
                for(int j=1;j<=n;j++){
                    tmp[j]+=f[0][j][i]*g[i];
                }
            }
        }
        memcpy(g,tmp,sizeof(tmp));
        for(int i=1;i<=n;i++)ans=max(ans,g[i]);
        printf("%.10lf\n",ans);
    }
}