• [Luogu4986] 逃离


    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;
    }
    
    
  • 相关阅读:
    Oracle EBS 更改物料说明后,在MTL_SYSTEM_ITEMS_B表中无变化
    Oracle EBS 复制用户职责
    Oracle EBS 多节点停应用
    Oracle EBS AR 更新客户
    Oracle EBS AP更新供应商地址
    Oracle EBS AP 供应商地点失效
    Oracle AP更新供应商
    string 从下标0 一直截到倒数第三位
    一个不错的spring 学习博客
    Filter及FilterChain的使用详解
  • 原文地址:https://www.cnblogs.com/YoungNeal/p/10300652.html
Copyright © 2020-2023  润新知