CF 963E Circles of Waiting

2023-05-16

- - - - 传送门
      - 题目大意
      - 思路
      - 参考代码

传送门

题目大意

在平面直角坐标系上，有一个神奇的点，一开始在 $(0, 0)$ 。每秒钟这个点都会随机移动：如果它在 $(x, y)$ ，下一秒它在 $(x - 1, y)$ 的概率是 $p_1$ ，在 $(x, y - 1)$ 的概率是 $p_2$ ，在 $(x + 1, y)$ 的概率是 $p_3$ ，在 $(x, y + 1)$ 的概率是 $p_4$ 。保证 $p_1 + p_2 + p_3 + p_4 = 1$ ，各次移动互不关联。

求出这个点移动至距离原点距离大于 $R \pod {R \le 50}$ 的点的期望步数。距离为欧几里得距离。结果对 $10^9 + 7$ 取模。

思路

唉，我太弱了，什么都不会，这道题这么显然需要解方程，但我还是不会做。

显然有：

fx,y=p1⋅fx−1,y+p2⋅fx,y−1+p3⋅fx+1,y+p4⋅fx,y+1+1 f x , y = p 1 ⋅ f x − 1 , y + p 2 ⋅ f x , y − 1 + p 3 ⋅ f x + 1 , y + p 4 ⋅ f x , y + 1 + 1

$f_{x, y} = p_1 \cdot f_{x - 1, y} + p_2 \cdot f_{x, y - 1} + p_3 \cdot f_{x + 1, y} + p_4 \cdot f_{x, y + 1} + 1$

总共有 $R^2$ 个点，因此我们立即得到一个 $O(R^6)$ 的算法（此处 $R$ 取 $100$ ），然后 T 掉。

感性地理解， $f_{x, y}$ 有系数的方程并不会太多，而且有系数的一定与它相邻。假设我们从上往下，从左往右编码，则有系数的项距自己的编码不会超过 $2r$ ，因此整个方程组对应的增广矩阵大概长这样：

对于每行，只需要在一个红框内消元即可：

红框的大小为 $2R$ ，因此时间复杂度为 $O(R^4)$ 。

参考代码

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <cassert>
#include <cctype>
#include <climits>
#include <ctime>
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <list>
#include <functional>
typedef long long LL;
typedef unsigned long long ULL;
using std::cin;
using std::cout;
using std::endl;
typedef int INT_PUT;
INT_PUT readIn()
{
    INT_PUT a = 0; bool positive = true;
    char ch = getchar();
    while (!(ch == '-' || std::isdigit(ch))) ch = getchar();
    if (ch == '-') { positive = false; ch = getchar(); }
    while (std::isdigit(ch)) { a = a * 10 - (ch - '0'); ch = getchar(); }
    return positive ? -a : a;
}
void printOut(INT_PUT x)
{
    char buffer[20]; int length = 0;
    if (x < 0) putchar('-'); else x = -x;
    do buffer[length++] = -(x % 10) + '0'; while (x /= 10);
    do putchar(buffer[--length]); while (length);
}

const int mod = int(1e9) + 7;
const int maxn = 55;
int r;
int p1, p2, p3, p4;
int sum;

inline int code(int x) { return x + 50; }
inline double dis(int x, int y) { return std::sqrt(std::pow(x, 2) + std::pow(y, 2)); }
std::vector<std::pair<int, int> > pt;
int idx[maxn * 2][maxn * 2];
#define IDX(x, y) idx[code(x)][code(y)]

LL power(LL x, int y)
{
    LL ret = 1;
    while (y)
    {
        if (y & 1) ret = ret * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ret;
}

int rect[8000][8000];

void Gauss()
{
    int d = r * 2;
    for (int i = 1; i <= pt.size(); i++)
    {
        LL inv = power(rect[i][i], mod - 2);
        for (int j = i + 1, to1 = 0; j <= pt.size() && to1 <= d; j++, to1++)
        {
            LL ratio = inv * rect[j][i] % mod;
            if (ratio)
            {
                for (int k = i, to2 = 0; k <= pt.size() && to2 <= d; k++, to2++)
                    rect[j][k] = ((LL)rect[j][k] - rect[i][k] * ratio) % mod;
                rect[j][0] = ((LL)rect[j][0] - rect[i][0] * ratio) % mod;
            }
        }
    }

    for (int i = pt.size(), to = IDX(0, 0); i >= to; i--)
    {
        LL inv = power(rect[i][i], mod - 2);
        int& t = rect[i][0];
        t = inv * t % mod;
        for (int j = i - 1; j >= to; j--)
            rect[j][0] = (rect[j][0] - (LL)rect[j][i] * t) % mod;
    }
}

void solve()
{
    for (int x = r; x >= -r; x--)
        for (int y = -r; y <= r; y++)
            if (dis(x, y) <= r)
            {
                pt.push_back(std::make_pair(x, y));
                IDX(x, y) = pt.size();
            }

    for (int x = r; x >= -r; x--)
        for (int y = -r; y <= r; y++)
            if (int f = IDX(x, y))
            {
                rect[f][f] = -1;
                if (int t = IDX(x - 1, y))
                    rect[f][t] = p1;
                if (int t = IDX(x, y - 1))
                    rect[f][t] = p2;
                if (int t = IDX(x + 1, y))
                    rect[f][t] = p3;
                if (int t = IDX(x, y + 1))
                    rect[f][t] = p4;
                rect[f][0] = -1;
            }

    Gauss();
    printOut((rect[IDX(0, 0)][0] + mod) % mod);
}

void run()
{
    r = readIn();
    p1 = readIn();
    p2 = readIn();
    p3 = readIn();
    p4 = readIn();
    sum = ((LL)p1 + p2 + p3 + p4) % mod;
    LL inv = power(sum, mod - 2);
    p1 = (LL)p1 * inv % mod;
    p2 = (LL)p2 * inv % mod;
    p3 = (LL)p3 * inv % mod;
    p4 = (LL)p4 * inv % mod;
    solve();
}

int main()
{
    run();
    return 0;
}

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)