Steven_Meng's Blog

前置芝士：
牛顿迭代
FFT

首先，我们将运动轨迹画出来：

发现$hdxrie$只能沿半径走，$Althen$只能沿竖直或横向走。

不妨将$Althen$的路径平移到两条垂直的线段上，总长不变。
设移动时间为$t$，由勾股定理，我们有$(A(x) \times t)^2=(B(x) \times t)^2+(C(x) \times t)^2$
即$A(x)^2=B(x)^2+C(x)^2$
发现$Len_a \le 10^5$，所以不能用朴素算法来求$A(x)^2$，考虑时间复杂度为$O(n \log n)$的快速傅里叶变换，我们可以很快地求出$A(x)^2$，$B(x)^2$和$C(x)^2$
最后，我们把$B(x)^2,C(x)^2$移项到左边，得到一个函数$A(x)^2-B(x)^2-C(x)^2=0$
我们可以用牛顿迭代求$A(x)^2-B(x)^2-C(x)^2=0$的所有解。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <vector>
#define MAXN 300005
#define MoHa 19260817
int tim=30;
const double eps=1e-10,PI=acos(-1.0);
using namespace std;
inline int iread(){
	int x=0,f=1;
	char ch=getchar();
	while (ch<'0'||ch>'9'){
		if (ch=='-') f=-1;
		ch=getchar();
	}
	while (ch<='9'&&ch>='0'){
		x=(x<<3)+(x<<1)+(ch^48);
		ch=getchar();
	}
	return x*f;
}
typedef vector<int> poly;
namespace Poly{
	struct complex{
		double x,y;
	};
	complex operator + (const complex &a,const complex &b){
		return (complex){a.x+b.x,a.y+b.y};
	}
	complex operator - (const complex &a,const complex &b){
		return (complex){a.x-b.x,a.y-b.y};
	}
	complex operator * (const complex &a,const complex &b){
		return (complex){a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};
	}
	complex cA[MAXN],cB[MAXN];
	int r[MAXN],pow2,n,m;
	inline void FFT(complex *A,int val,int len){
		for (int i=0;i<len;++i){
			if (i<r[i]) swap(A[i],A[r[i]]);
		}
	    for (int i=1;i<len;i<<=1){
	    	complex Wn=(complex){cos(PI/i),val*sin(PI/i)};
	    	for (int j=0;j<len;j+=(i<<1)){
	    		complex t=(complex){1,0};
	    		for (int k=0;k<i;++k,t=t*Wn){
	    			complex x=A[j+k],y=t*A[i+j+k];
	    			A[j+k]=x+y;
	    			A[i+j+k]=x-y;
				}
			}
		}
		if (val==-1) for (int i=0;i<len;++i) A[i].x/=len;
	}
	inline void init_fft(poly a){
		for (int i=0;i<a.size();++i){
			cA[i].x=cB[i].x=(double)a[i];
			cA[i].y=cB[i].y=0;
		}
	}
	inline poly mul(poly a){
		init_fft(a);
		int n=a.size()-1,m=a.size()-1;
		a.resize(m+2);
		int L=0;
		for (m+=n,n=1;n<=m;n<<=1) L++;
		for (int i=0;i<n;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(L-1));
		FFT(cA,1,n),FFT(cB,1,n);
		for (int i=0;i<=n;++i) cA[i]=cA[i]*cB[i];
		FFT(cA,-1,n);
		for (int i=0;i<a.size();++i) a[i]=(int)(cA[i].x+0.1);
		return a;
	}
}
using namespace Poly;
poly SA,SB,SC,C1;
inline void read(poly &A,int len){
	for (int i=0;i<=len;++i) A.push_back(iread());
}
inline double F(poly &R,double x){
	double ans=0;
	for (int i=R.size()-1;i>=0;--i) ans=ans*x+(double)R[i];
	return ans;
}
inline void qd(poly &F){
	C1.resize(F.size()-1);
	for (int i=0;i<C1.size();++i){
		C1[i]=F[i+1]*(i+1);
	}
}
double L,R;
inline double Newton(double x){
	qd(SC);
	double c;
	while (true){
		--tim;
		c=F(SC,x);
		if (fabs(c)<eps) return x;
		x=x-c/F(C1,x);
		x=max(x,L),x=min(x,R);
		if (!tim) return 0;
	}
}
int main(){
	int la,lb,lc;
	scanf("%d%d%d%lf%lf",&la,&lb,&lc,&L,&R);
	read(SA,la),read(SB,lb),read(SC,lc);
	SA=mul(SA),SB=mul(SB),SC=mul(SC);
	SC.resize(max(max(SA.size(),SB.size()),SC.size()));
	for (int i=0;i<SC.size();++i){
		SC[i]-=SA[i]+SB[i];
	}
	double ans=Newton((L+R)/2.00);
	if (!tim) printf("Inconsistent!\n");
	else printf("%.8lf\n",ans);
}