LJJ的电阻?格

题目概述

在这里插入图片描述

题解

~~什么**出题人，在这里出物理题，你当人人都是物竞的是吧。~~

首先，我们可以一个一个点地来想，不过我们首先可以将所有的点分成两大类 $x = y$ 与 $x\not = y$ 两大类。
对于 $x = y$ 的部分，观察到所有的 $x, y$ 范围都是远大于 $x\not = y$ 部分的，而它又给了我们一个公式，显然是要让我们通过这个公式算它的答案。
对于 $x = 100000, y = 100000$ 的点，显然可以按它的公式，算出 $\sum_{i=1}^x\frac{1}{2i-1}$ ，再乘上 $\frac{2}{\pi}$ 就行了。
而 $x,y\in[10^{10},10^{15}]$ 的部分显然是让我们通过更快的方式去计算这个值。
$\lim_{x\rightarrow +\infty}\sum_{i=1}^x\frac{1}{2i-1}=\lim_{x\rightarrow +\infty}\sum_{i=1}^{2x}\frac{1}{i}-\frac{1}{2}\lim_{x\rightarrow +\infty}\sum_{i=1}^x\frac{1}{i}$
而我们又知道 $\lim_{x\rightarrow \infty}\sum_{i=1}^{x}\frac{1}{i}$ 是与 $\ln n$ 比较接近的，由莱昂哈德·欧拉告诉我们它们的差值是趋近于 $\gamma=0.577215664901532860606512090082402431042159335......$ ，于是我们的答案就是 $\ln(2x)-\frac{1}{2}\ln (x)+\frac{\gamma}{2}$ 。

而对于 $x\not = y$ 的部分，我们观察到 $x, y$ 都相当的小的，我们可以用稍微暴力的方法求解。
对于 $x = 1, y = 0$ 的部分，学过高中物理的人应该比较容易算出答案。
我们要计算的是 $(0, 0) ? > (0, 1)$ 这一条边的等效电阻，可以通过计算该边上的电流计算出。
该电流，我们可以通过施加电场构造出来。
施加第一个电场，使得有 $I$ 的电流从点 $a$ 进入，向无限蔓延出去，显然，这样的话我们的电场应该达到如下的效果，
在这里插入图片描述
显然，通过对称性，点 $a$ 向四个方向中每个方向流去的电流都应该是 $\frac{I}{4}$ ，施加在 $a ? > b$ 的边上的电流应该是 $\frac{I}{4}$ 。
我们再在点 $b$ 上施加一个类似但反向的电场，使得电流从四面八方涌来，从 $b$ 流出一个 $I$ 的电流，根据电场的可叠加性，这样的话，整张图就不会有电荷的堆积，而是形成了一个大小为 $I$ 的回路。
同样，这个电场也会在 $a ? > b$ 的边上有 $\frac{I}{4}$ 的电流。这条边总的电流就是 $\frac{I}{2}$ ，那么我们总电动势为 $\frac{UI}{2}$ ，等效电阻为 $\frac{R}{2}$ 。

我们记 $I (x, y)$ 表示从 $(0, 0)$ 到点 $(x, y)$ 的电流大小，我们已经轻易地求出了 $I (1, 0)$ 与 $I (x, x)$ ，我们要尝试通过这些算出我们剩下的 $I (x, y)$ 。
显然，从一个点流入的电流与从其流出的电流应该大小是相同的，否则它就要累积电荷了。
而根据对称性，每个点向四周流去的电流应该是相同的，那么我们就有 $4 I (x, y) = I (x ? 1, y) + I (x, y ? 1) + I (x + 1, y) + I (x, y + 1)$ 成立。
我们可以通过我们的已有信息算出未知信息。
在这里插入图片描述
我们可以不断枚举 $x + y = k$ 上的点，一层一层下去，按照向量的方向构造出下一层的向量。
根据上面的式子，在我们已知 $I (x, y), I (x ? 1, y), I (x, y ? 1)$ 时，再知道 $I (x + 1, y)$ 与 $I (x, y + 1)$ 中的一者就可以知道另一者。
对于中间 $y = x$ 上的点，其两侧的 $I (x + 1, y)$ 与 $I (x, y + 1)$ 显然是等效的，所以我们只需要知道前三者就行。 $I (x, y)$ 我们已预先知道， $I (x ? 1, y)$ 与 $I (x, y ? 1)$ 也已在上一层得到，我们显然可以轻松的知道 $I (x + 1, y)$ 与 $I (x, y + 1)$ 。
而往该中线的外侧扩张时，对于下一个点，我们也已经构造出两者中的一者，自然也能得到另一个。
到了坐标轴上，我们可以将坐标轴这边的翻转过去得到另一侧的值，自然可以构造出顶上的点。
按这个顺序就可以得到我们的 $I (x, y)$ 了。
于是乎，我们就轻松地解决该问题了。

时间复杂度 $O([x\not =y]\max(x,y)^2)$ 。
不过按题解说法， $x, y$ 大了就不精确了，所以只开到了 $15$ 。

源码

#include<bits/stdc++.h>
using namespace std;
#define MAXN 80005
#define lowbit(x) (x&-x)
#define reg register
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
#define lson (rt<<1)
#define rson (rt<<1|1)
typedef long long LL;
typedef unsigned long long uLL; 
typedef long double ld;
typedef pair<int,int> pii;
const int INF=0x3f3f3f3f;
const int mo=998244353;
const int mod=1e5+3;
const int inv2=499122177;
const int jzm=2333;
const int zero=2000;
const int n1=2000;
const int orG=3,ivG=332748118;
const long double Pi=acos(-1.0);
const double eps=1e-12;
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
template<typename _T>
void read(_T &x){
	_T f=1;x=0;char s=getchar();
	while(s>'9'||s<'0'){if(s=='-')f=-1;s=getchar();}
	while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
	x*=f;
}
template<typename _T>
void print(_T x){if(x<0){x=(~x)+1;putchar('-');}if(x>9)print(x/10);putchar(x%10+'0');}
int gcd(int a,int b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1)t=1ll*t*a%p;a=1ll*a*a%p;s>>=1;}return t;}
LL x,y;double ans,dp[50][50];
signed main(){
	freopen("grid.in","r",stdin);
	freopen("grid.out","w",stdout);
	read(x);read(y);
	if(x==1&&y==0){puts("0.500000");return 0;}
	if(x==y){
		if(x<=100000)for(int i=1;i<=x;i++)ans+=1.0/(i+i-1);
		else ans=log(2.0*x)-0.5*(log(1.0*x)-0.57721566490153286060);
		ans*=2.0/Pi;
		printf("%.6f\n",ans);
		return 0;
	}
	dp[0][0]=0;dp[1][0]=dp[0][1]=0.5;
	for(int i=1;i<=max(x,y);i++)dp[i][i]=dp[i-1][i-1]+2.0/Pi/(i+i-1);
	for(int k=1;k<=x+y;k++)
		for(int ix=k/2,iy=(k+1)/2,jx=(k+1)/2,jy=k/2;ix>=0&&jy>=0;ix--,iy++,jx++,jy--){
			if(ix==iy)dp[ix+1][iy]=dp[ix][iy+1]=2.0*dp[ix][iy]-dp[ix-1][iy];
			else if(ix)dp[ix][iy+1]=4.0*dp[ix][iy]-dp[ix-1][iy]-dp[ix][iy-1]-dp[ix+1][iy];
			else dp[ix][iy+1]=4.0*dp[ix][iy]-dp[ix][iy-1]-2.0*dp[ix+1][iy];
			if(jx==jy)dp[jx+1][jy]=dp[jx][jy+1]=2.0*dp[jx][jy]-dp[jx-1][jy];
			else if(jy)dp[jx+1][jy]=4.0*dp[jx][jy]-dp[jx-1][jy]-dp[jx][jy-1]-dp[jx][jy+1];
			else dp[jx+1][jy]=4.0*dp[jx][jy]-dp[jx-1][jy]-2.0*dp[jx][jy+1];
		}
	printf("%.6f\n",dp[x][y]);
	return 0;
}