[Luogu4986] 逃离

Description

给定次数为 \(n\) 的函数 \(A(x),B(x),C(x)\)，求 \(A^2(x)+B^2(x)-C^2(x)\) 在 \([L,R]\) 的零点。\(n\leq 10^5,-3\le L,R\le 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;
}