Steven_Mengの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$的所有解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#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);
}

 评论