Description
给定次数为 (n) 的函数 (A(x),B(x),C(x)),求 (A^2(x)+B^2(x)-C^2(x)) 在 ([L,R]) 的零点。(nleq 10^5,-3le L,Rle 3)。
Solution
先对这三个函数 (FFT) 一下,然后乘起来求出新函数 (f(x)),牛顿迭代一下就没了。
为了防止以后忘记牛顿迭代怎么写,在这里mark一下。没错我就是马来人
牛顿迭代本质上就是对于原函数 (f(x)) 构造一个新函数 (g(x)) 使得在 (x_i) 这个点的时候,(f(x_i)=g(x_i)land f'(x_i)=g'(x_i))。也就是让原函数和一阶导都相等。构造出来 (g) 之后每次让 (y=0) 然后取 (x) 的新解。这样迭代几次就会发现 (f(x)) 趋向于 (0) 了。
哦对公式就是 (x_n=x_{n-1}-frac{f(x_{n-1})}{f'(x_{n-1})})
Code
#include<bits/stdc++.h>
using std::min;
using std::max;
using std::swap;
using std::vector;
typedef double db;
typedef long long ll;
#define pb(A) push_back(A)
#define pii std::pair<int,int>
#define all(A) A.begin(),A.end()
#define mp(A,B) std::make_pair(A,B)
const int N=4e5+5;
const db Pi=acos(-1);
db L,R;
int la,lb,lc,n;
int lim,rev[N];
db a[N],b[N],c[N];
db a0[N],b0[N],c0[N];
struct cp{
db x,y;
cp(db a=0,db b=0){x=a,y=b;}
friend cp operator+(cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
friend cp operator-(cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
friend cp operator*(cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
}f[N],g[N];
void fft(cp *f,int opt){
for(int i=1;i<lim;i++) if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1){
cp tmp=cp(cos(Pi/mid),sin(Pi/mid)*opt);
for(int R=mid<<1,j=0;j<lim;j+=R){
cp w=cp(1,0);
for(int k=0;k<mid;k++,w=w*tmp){
cp x=f[j+k],y=f[j+k+mid]*w;
f[j+k]=x+y;f[j+k+mid]=x-y;
}
}
} if(opt<0)
for(int i=0;i<lim;i++) f[i].x/=lim;
}
void merge(db *a,db *b,int n){
lim=1;while(lim<=n+n) lim<<=1;
for(int i=1;i<lim;i++) rev[i]=(rev[i>>1]>>1)|(i&1?lim>>1:0);
for(int i=0;i<=n;i++) f[i]=cp(a[i],0);
for(int i=n+1;i<lim;i++) f[i]=cp(0,0);
fft(f,1);
for(int i=0;i<lim;i++) f[i]=f[i]*f[i];
fft(f,-1);
for(int i=0;i<lim;i++) b[i]=f[i].x;
}
int getint(){
int X=0,w=0;char ch=getchar();
while(!isdigit(ch))w|=ch=='-',ch=getchar();
while( isdigit(ch))X=X*10+ch-48,ch=getchar();
if(w) return -X;return X;
}
db ds(db x){
db now=1,ans=0;
for(int i=1;i<=n;i++)
ans+=a[i]*now*i,now*=x;
return ans;
}
db yh(db x){
db now=1,ans=0;
for(int i=0;i<=n;i++)
ans+=a[i]*now,now*=x;
return ans;
}
void solve(db x){
for(int i=1;i<=30;i++){
db now=yh(x);
if(fabs(now)<1e-9) return printf("%.8lf",x),void();
x=x-now/ds(x);
x=max(x,L);x=min(x,R);
} puts("Inconsistent!");
}
signed main(){
la=getint(),lb=getint(),lc=getint();scanf("%lf%lf",&L,&R);
for(int i=0;i<=la;i++) scanf("%lf",&a0[i]);
for(int i=0;i<=lb;i++) scanf("%lf",&b0[i]);
for(int i=0;i<=lc;i++) scanf("%lf",&c0[i]);
merge(a0,a,la),merge(b0,b,lb),merge(c0,c,lc);
n=max(la,max(lb,lc))<<1;
for(int i=0;i<=n;i++) a[i]+=b[i]-c[i];
solve((L+R)/2); return 0;
}