求助!UKE怎么解?
查看原帖
求助!UKE怎么解?
479909
708zz楼主2023/7/23 14:50
#include <bits/stdc++.h>
using namespace std;

inline int read(){
	int sum=0,f=0;
	char ch=getchar();
	for(;!isdigit(ch);ch=getchar()){
		f |= (ch=='-');
	}
	for(;isdigit(ch);ch=getchar()){
		sum = ((sum<<3) + (sum<<1) + (ch^48));
	}
	return f?-sum:sum;
}

typedef long long ll;

const int maxn = 100;
const int mod = 1e9+7;

ll n,m,r,a1,a2,a3,a4,tot,an[maxn],mp[maxn][maxn],ans;
bool flag[maxn*maxn],ed[maxn*maxn];

struct vec{
	ll x,y;
	friend bool operator < (vec x,vec y) {
		return x.x > y.x;
	}
};

priority_queue <vec> q;

struct node{
	ll a[maxn];
	friend node operator * (ll x,node y){
		for(int i=0;i<=2*r+1;i++){
			y.a[i] = y.a[i]*x%mod;
		}
		return y;
	}
	friend node operator + (node y,ll x){
		y.a[2*r+1] = (y.a[2*r+1]+x)%mod;
		return y;
	}
	friend node operator + (node x,node y){
		for(int i=0;i<=2*r+1;i++){
			x.a[i] = (x.a[i]+y.a[i])%mod;
		}
		return x;
	}
}f[maxn*maxn];

ll power(ll x,ll y){
	ll res = 1;
	while(y){
		if(y&1){
			res = res*x%mod;
		}
		x = x*x%mod;
		y >>= 1;
	}
	return res;
}

ll inv(ll x){
	return power(x,mod-2);
}

ll mark(ll x,ll y){
	return (y+r+1)*n+(x+r+2);
}

void solve(){
	for(int i=0;i<=2*r;i++){
		if(!mp[i][i]){
			ll kt = 0;
			for(int j=i+1;j<=2*r;j++){
				if(mp[j][i]){
					for(int k=i;k<=2*r+1;k++){
						swap(mp[i][k] , mp[j][k]);
					}
					kt = 1;
					break;
				}
			}
			if(!kt){
				return ;
			}
		}
		ll in = inv(mp[i][i]);
		for(int j=i+1;j<=2*r;j++){
			ll op = (mod-mp[j][i])*in%mod;
			for(int k=i;k<=2*r+1;k++){
				mp[j][k] = (mp[j][k] + op*mp[i][k]%mod+mod)%mod;
			}
		}
	}
	for(int i=2*r;i>=0;i--){
		an[i] = mp[i][2*r+1]*inv(mp[i][i])%mod;
		for(int j=0;j<=i-1;j++){
			mp[j][2*r+1] = (mp[j][2*r+1]-an[i]*mp[j][i]%mod+mod)%mod;
		}
	}
}

int main(){
	r=read();
	a1=read(); a2=read(); a3=read(); a4=read();
	n = 2*r+3;
	ll pg = inv(a1+a2+a3+a4);
	a1 = a1*pg%mod;
	a2 = a2*pg%mod;
	a3 = a3*pg%mod;
	a4 = a4*pg%mod;
	for(int i=-r;i<=r;i++){
		for(int j=-r;j<=r;j++){
			if(i*i+j*j <= r*r){
				flag[mark(i,j)] = 1;
			}
		}
	}
	for(int j=-r;j<=r;j++){
		for(int i=-r;i<=0;i++){
			if(flag[mark(i,j)]){
				f[mark(i,j)].a[j+r] = 1;
				q.push(vec{i,j});
				break;
			}
		}
	}
	for(int j=-r;j<=r;j++){
		for(int i=r;i>=-r;i--){
			if(flag[mark(i,j)]){
				ed[mark(i+1,j)] = 1;
				break;
			}
		}
	}
	while(!q.empty()){
		ll X = q.top().x, Y = q.top().y;
		q.pop();
		if(!flag[mark(X+1,Y)] && !ed[mark(X+1,Y)]){
			continue;
		}
		f[mark(X+1,Y)] = inv(a3)*((f[mark(X,Y)]+(mod-1)) + ((mod-a1)*f[mark(X-1,Y)]) + ((mod-a2)*f[mark(X,Y-1)]) + ((mod-a4)*f[mark(X,Y+1)]));
		if(ed[mark(X+1,Y)]){
			for(int i=0;i<=2*r+1;i++){
				mp[Y+r][i] = f[mark(X+1,Y)].a[i];
			}
			mp[Y+r][2*r+1] = mod-mp[Y+r][2*r+1];
		}else {
			q.push(vec{X+1,Y});
		}
	}
	solve();
	an[2*r+1] = 1;
	for(int i=0;i<=2*r+1;i++){
		ans = (ans + f[mark(0,0)].a[i]*an[i]%mod)%mod;
	}
	printf("%lld\n",ans);
	return 0;
}
2023/7/23 14:50
加载中...