Codeforces989E. A Trance of Nightfall

$n \leq 200$个平面上的点，$q \leq 200$次询问：重复操作$m \leq 10000$次，到达点$x$的概率最大是多少。操作：一开始选点$P$，不一定要是给定点，可以是平面上任一点。然后，选一条穿过给定点至少两个点且穿过$P$的直线$l$，若有多条，等概率选一条；选中一条后，把$P$点移动到这条直线上的某个初始给定点，若有多个等概率选。

以为是不可做几何题，其实挺好玩的

如果知道两个点之间一次性直接到达的概率$A(x,y)$，那么记$f(i,j)$为走$i$步到点$j$的概率，$f(i,j)=\sum_{k=1}^{n}A(i,k)f(i-1,k)$，是个矩阵乘法，可以快速幂优化。后面$m-1$步可以用矩阵算，那第一步应该怎么选初始点呢？可以发现，固定了终点，用矩阵乘法推出了$f(m,i)$后，每条直线的贡献是固定的，而如果我们选的$P$点在几个直线的交汇处，得到的便是这几条直线贡献的平均值，这不如直接选一条贡献最大的直线（这里的选是说，把点P放在这条线上），而这总是可以做到的。因此输出答案最大的一条直线的答案即可。

有个问题，如果每次询问都要快速幂一次，那不久$qn^3logn$了吗？别担心，$A^{2^i}$是可以预处理的，等会要做乘法的时候只需拿$f$和某些$A^{2^i}$乘就好了，而一个一维数组$f$和矩阵的乘法是$n^2$的，如此少一个$n$，可通过。

通过个鬼。直线运算精度误差坑了我一个半小时。。这里上一个标程无误差版的板子：

struct line {
    // ax + by + c = 0
    int a, b, c;

    line() : a(0), b(0), c(0) { }
    line(const point &p, const point &q)
        : a(p.y - q.y)
        , b(q.x - p.x)
        , c(q.y * p.x - q.x * p.y)
    {
        int g = gcd(gcd(a, b), c);
        if (g != 1) { a /= g; b /= g; c /= g; }
        if (a < 0) { a = -a; b = -b; c = -c; }
    }

    inline bool contains(const point &p) {
        return a * p.x + b * p.y + c == 0;
    }

    // For sorting & duplicate removing
    inline bool operator < (const line &other) const {
        return a != other.a ? a < other.a :
            b != other.b ? b < other.b : c < other.c;
    }
    inline bool operator == (const line &other) const {
        return a == other.a && b == other.b && c == other.c;
    }
};

以及我改longdouble终于通过的辣鸡版本：

#include<iostream>
#include<cstring>
#include<cstdio>
#include<math.h>
//#include<complex>
//#include<set>
//#include<queue>
#include<vector>
#include<algorithm>
#include<stdlib.h>
using namespace std;

#define double long double
#define LL long long
int qread()
{
    char c; int s=0,f=1; while ((c=getchar())<'0' || c>'9') (c=='-') && (f=-1);
    do s=s*10+c-'0'; while ((c=getchar())>='0' && c<='9'); return s*f;
}

//Pay attention to '-' , LL and double of qread!!!!

int n;
#define maxn 40011
#define maxm 40011
#define eps 1e-10
bool equ(double a,double b) {return fabs(a-b)<eps;}
bool neq(double a,double b) {return fabs(a-b)>eps;}
bool le(double a,double b) {return a-b<-eps;}
bool leq(double a,double b) {return a-b<eps;}
bool ge(double a,double b) {return a-b>eps;}
bool geq(double a,double b) {return a-b>-eps;}

struct Point{int x,y;}p[maxn];
struct Line
{
    double k,b;
    void make(const Point &A,const Point &B)
    {
        if (A.x==B.x) {k=1e18; b=A.x; return;}
        k=1.0*(A.y-B.y)/(A.x-B.x);
        b=A.y-k*A.x;
    }
    bool operator < (const Line &l)
    {
        if (neq(l.k,k)) return le(k,l.k);
        return le(b,l.b);
    }
    bool operator == (const Line &l) {return equ(l.k,k) && equ(l.b,b);}
    bool operator != (const Line &l) {return !((*this)==l);}
}l[maxm]; int len=0;
vector<int> vl[maxm],vp[maxn];

struct Mat
{
    double a[211][211];
    Mat() {memset(a,0,sizeof(a));}
    Mat operator * (const Mat &b)
    {
        Mat ans;
        for (int i=1;i<=n;i++)
            for (int j=1;j<=n;j++)
                for (int k=1;k<=n;k++)
                    ans.a[i][j]+=a[i][k]*b.a[k][j];
        return ans;
    }
}base[20];

double f[maxn],g[maxn];
void mul(double *f,Mat b)
{
    memset(g,0,sizeof(g));
    for (int j=1;j<=n;j++)
        for (int k=1;k<=n;k++)
            g[j]+=f[k]*b.a[j][k];
    for (int j=1;j<=n;j++) f[j]=g[j];
}

int main()
{
    n=qread();
    for (int i=1;i<=n;i++) {p[i].x=qread(); p[i].y=qread();}
    for (int i=1;i<=n;i++)
        for (int j=i+1;j<=n;j++)
        {
            len++;
            l[len].make(p[i],p[j]);
        }
    sort(l+1,l+1+len);
    l[0].k=l[len+1].k=-1e18; l[0].b=l[len+1].b=-1e18;
    {
        int j=0;
        for (int i=1;i<=len;i++) if (l[i]!=l[i-1]) l[++j]=l[i];
        len=j;
    }
    
    for (int i=1;i<=n;i++)
        for (int j=1;j<=len;j++)
            if ((neq(l[j].k,1e18) && equ(l[j].k*p[i].x+l[j].b,p[i].y))
            || (equ(l[j].k,1e18) && equ(p[i].x,l[j].b)))
                vl[j].push_back(i),vp[i].push_back(j);
    
    for (int i=1;i<=n;i++)
        for (int j=0,to=vp[i].size();j<to;j++)
        {
            int u=vp[i][j];
            for (int k=0,too=vl[u].size();k<too;k++)
            {
                int x=vl[u][k];
                base[0].a[i][x]+=1.0/(to*too);
            }
        }
    
    for (int i=1;i<=17;i++) base[i]=base[i-1]*base[i-1];
    int lq=qread(),x,y;
    while (lq--)
    {
        x=qread(); y=qread(); y--;
        memset(f,0,sizeof(f)); f[x]=1;
        for (int i=0;i<=17;i++) if ((y>>i)&1) mul(f,base[i]);
        double ans=0;
        for (int i=1;i<=len;i++)
        {
            double tmp=0;
            for (int j=0,to=vl[i].size();j<to;j++) tmp+=f[vl[i][j]];
            ans=max(ans,tmp/vl[i].size());
        }
        cout<<ans<<endl;
    }
    return 0;
}