多项式快速插值

200+行的多项式板子题真爽啊

给定\(n\)个点的点值\((x_i,y_i)\),求这\(n\)个点确定的\(n-1\)次多项式

\(n\le 10^5\)

前置知识:

多项式多点求值

拉格朗日插值

微积分基础

首先我们有一个\(n^2\)的拉格朗日插值法

\[f(x)=\sum\limits_{i=1}^{n}y_i\prod\limits_{i\ne j}\frac{x-x_j}{x_i-x_j}\]

然后我们学习一个WC2017挑战就过了

考虑优化,我们知道这个形式它很死,把它变成重心插值

\[f(x)=\sum\limits_{i=1}^{n}\frac{y_i}{\prod\limits_{i\ne j}x_i-x_j}\prod\limits_{i\ne j}(x-x_j)\]

发现除了一个常数\(y_i\)剩下的就是\(\prod\limits_{i\ne j}(x-x_i)\)的形式了,我们设它为\(g(x)\)

那么前面那一项的分母可以表示为\(\frac{g(x)}{x-x_i}(x=x_i)\)

发现\(x=x_i\)时分子分母都是\(0\)

根据洛必达法则

\[\lim\limits_{x→x_i}\frac{g(x_i)}{x-x_i}=\lim\limits_{x→x_i}\frac{g'(x_i)}{(x-x_i)'}=\lim\limits_{x→x_i}g'(x)\]

分治算出\(g(x)\)后多点求值算出\(g'(x_i)\)

然后拿\(y_i\)除一下,前一项就搞定了

\[f_{l,r}=\sum\limits_{i=l}^{r}\frac{y_i}{g'(x_i)}\prod\limits_{j=l,j\ne i}^{r}(x-x_j)\]

\[=\prod\limits_{j=mid+1}^{r}(x-x_j)\sum\limits_{i=l}^{mid}\frac{y_i}{g'(x_i)}\prod\limits_{j=l,j\ne i}^{mid}(x-x_j)+\prod\limits_{j=l}^{mid}(x-x_j)\sum\limits_{i=mid+1}^{r}\frac{y_i}{g'(x_i)}\prod\limits_{j=mid+1,j\ne i}^{r}(x-x_j)\]

\[=\prod\limits_{i=mid+1}^{r}(x-x_i)f_{l,mid}+\prod\limits_{i=l}^{mid}(x-x_i)f_{mid+1,r}\]

\[=g_{mid+1,r}f_{l,mid}+g_{l,mid}f_{mid+1,r}\]

分治求解

这板子让我再打一遍都打不出来

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define ls(p) (p<<1)
#define rs(p) (p<<1|1)
#define eps (1e-8)
    inline int read()
    {
        int x=0;char ch,f=1;
        for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
        if(ch=='-') f=0,ch=getchar();
        while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
        return f?x:-x;
    }
    const int N=1e5+10,M=18,mod=998244353;
    int n,m,tot,limit,len;
    int xx[N],yy[N],val[N];
    int a[N],b[N],c[N],rr[N],cc[N<<2];
    int ra[N],rb[N<<2],irb[N<<2];
    int f[N*M<<1],stf[N<<2],enf[N<<2];
    int g[N*M<<1],stg[N<<2],eng[N<<2];
    int h[N],sth[N],enh[N],f3[N];
    int f1[N<<2],f2[N<<2],w[21][N<<2],pos[N<<2];
    inline int fast(int x,int k)
    {
        int ret=1;
        while(k)
        {
            if(k&1) ret=ret*x%mod;
            x=x*x%mod;
            k>>=1;
        }
        return ret;
    }
    inline int add(int x,const int &y)
    {
        x+=y;
        return x>mod?x-mod:x;
    }
    inline int del(int x,const int &y)
    {
        x-=y;
        return x<0?x+mod:x;
    }
    inline void init(int x)
    {
        limit=1,len=0;
        while(limit<x) limit<<=1,++len;
        for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
    }
    inline void ntt(int limit,int *a,int inv)
    {
        for(int i=0;i<limit;++i)
            if(i<pos[i]) swap(a[i],a[pos[i]]);
        for(int mid=1,t=1;mid<limit;mid<<=1,++t)
        {
            for(int r=mid<<1,j=0;j<limit;j+=r)
            {
                for(int k=0;k<mid;++k)
                {
                    int x=a[j+k],y=w[t][k]*a[j+k+mid]%mod;
                    a[j+k]=add(x,y);
                    a[j+k+mid]=del(x,y);
                }
            }
        }
        if(inv) return;
        inv=fast(limit,mod-2);reverse(a+1,a+limit);
        for(int i=0;i<limit;++i) a[i]=a[i]*inv%mod;
    }
    inline void NTT(int *a,int *b,int limit)
    {
        ntt(limit,a,1);ntt(limit,b,1);
        for(int i=0;i<limit;++i) a[i]=a[i]*b[i]%mod;
        ntt(limit,a,0);
    }
    inline void poly_inv(int limit,int len,int *a,int *b)
    {
        if(limit==1){b[0]=fast(a[0],mod-2);return;}
        poly_inv(limit>>1,len-1,a,b);
        for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
        for(int i=0;i<(limit>>1);++i) cc[i]=a[i];
        for(int i=limit>>1;i<limit;++i) cc[i]=0;
        ntt(limit,cc,1);ntt(limit,b,1);
        for(int i=0;i<limit;++i) b[i]=((2-cc[i]*b[i]%mod)+mod)%mod*b[i]%mod;
        ntt(limit,b,0);
        for(int i=limit>>1;i<limit;++i) b[i]=0;
    }
    inline void make(int l,int r,int p)
    {
        if(l==r)
        {
            g[stg[p]=++tot]=mod-xx[l];
            g[eng[p]=++tot]=1;
            return;
        }
        int mid=(l+r)>>1;
        make(l,mid,ls(p));
        make(mid+1,r,rs(p));
        int na=eng[ls(p)]-stg[ls(p)]+1;
        int nb=eng[rs(p)]-stg[rs(p)]+1;
        init(na+nb);
        for(int i=0;i<na;i++) f1[i]=g[stg[ls(p)]+i];
        for(int i=na;i<limit;i++) f1[i]=0;
        for(int i=0;i<nb;i++) f2[i]=g[stg[rs(p)]+i];
        for(int i=nb;i<limit;i++) f2[i]=0;
        NTT(f1,f2,limit);
        stg[p]=tot+1;
        na+=nb-1;
        for(int i=0;i<na;++i) g[++tot]=f1[i];
        eng[p]=tot;
    }
    inline void solve(int l,int r,int p,int fa)
    {
        int na=enf[fa]-stf[fa],nb=eng[p]-stg[p];
        if(na>=nb)
        {
            int nc=na-nb;
            for(int i=0;i<=na;++i) a[i]=f[stf[fa]+i];
            for(int i=0;i<=nb;i++) b[i]=g[stg[p]+i];
            for(int i=0;i<=nc;i++) ra[i]=a[na-i];
            for(int i=0;i<=nb;i++) rb[i]=b[nb-i];
            for(int i=nc+1;i<=nb;i++) rb[i]=0;
            init(nc*2+2);
            for(int i=nb+1;i<limit;i++) rb[i]=0;
            for(int i=0;i<limit;i++) irb[i]=0,f1[i]=0;
            poly_inv(limit,len,rb,irb);
            for(int i=0;i<=nc;i++) f1[i]=ra[i],f2[i]=irb[i];
            for(int i=nc+1;i<limit;i++) f1[i]=f2[i]=0;
            NTT(f1,f2,limit);
            for(int i=0;i<=nc;i++) c[nc-i]=f1[i];
            for(int i=nc+1;i<nb;i++) c[i]=0;
            init(nb<<1);
            for(int i=0;i<nb;i++) f1[i]=b[i],f2[i]=c[i];
            for(int i=nb;i<limit;i++) f1[i]=0,f2[i]=0;
            NTT(f1,f2,limit);
            for(int i=0;i<nb;i++) rr[i]=(a[i]-f1[i]+mod)%mod;
            while(nb>1 && !rr[nb-1]) nb--;
            stf[p]=tot+1;
            for(int i=0;i<nb;i++) f[++tot]=rr[i];
            enf[p]=tot;
        }
        else
        {
            stf[p]=tot+1;
            for(int i=stf[fa];i<=enf[fa];++i) f[++tot]=f[i];
            enf[p]=tot;
        }
        if(l==r)
        {
            val[l]=f[stf[p]];
            return;
        }
        int mid=(l+r)>>1;
        solve(l,mid,ls(p),p);
        solve(mid+1,r,rs(p),p);
    }
    inline void work(int l,int r,int p)
    {
        if(l==r) return;
        int mid=(l+r)>>1;
        work(l,mid,ls(p));
        work(mid+1,r,rs(p));
        int na=enh[l]-sth[l]+1,nb=eng[rs(p)]-stg[rs(p)]+1;

        init(na+nb);
        for(int i=0;i<na;i++) f1[i]=h[sth[l]+i];
        for(int i=na;i<limit;i++) f1[i]=0;
        for(int i=0;i<nb;i++) f2[i]=g[stg[rs(p)]+i];
        for(int i=nb;i<limit;i++) f2[i]=0;
        NTT(f1,f2,limit);

        na+=nb-1;
        for(int i=0;i<na;i++) f3[i]=f1[i];
        for(int i=na;i<limit;i++) f3[i]=0;

        na=enh[mid+1]-sth[mid+1]+1,nb=eng[ls(p)]-stg[ls(p)]+1;
        for(int i=0;i<na;i++) f1[i]=h[sth[mid+1]+i];
        for(int i=na;i<limit;i++) f1[i]=0;
        for(int i=0;i<nb;i++) f2[i]=g[stg[ls(p)]+i];
        for(int i=nb;i<limit;i++) f2[i]=0;
        NTT(f1,f2,limit);
        na+=nb-1;
        for(int i=0;i<na;i++) h[sth[l]+i]=f3[i]+f1[i]>=mod?f3[i]+f1[i]-mod:f3[i]+f1[i];
        enh[l]=sth[l]+na-1;
    }
    inline void main()
    {
        n=read();
        for(int mid=1,t=1;mid<400000;mid<<=1,++t)
        {
            w[t][0]=1;int Wn=fast(3,(mod-1)/(mid<<1));
            for(int k=1;k<mid;++k)
            {
                w[t][k]=w[t][k-1]*Wn%mod;
            }
        }
        for(int i=1;i<=n;++i) xx[i]=read(),yy[i]=read();
        make(1,n,1);
        m=eng[1]-stg[1];
        for(int i=0;i<=m;++i) f1[i]=g[stg[1]+i];
        for(int i=0;i<m;++i) f1[i]=f1[i+1]*(i+1)%mod;
        f1[m--]=0;
        stf[tot=0]=1;
        for(int i=0;i<=m;++i) f[enf[0]=++tot]=f1[i];
        solve(1,n,1,0);
        for(int i=1;i<=n;i++) val[i]=yy[i]*fast(val[i],mod-2)%mod;
        tot=0;
        for(int i=1;i<=n;i++)
        {
            sth[i]=enh[i]=++tot;
            h[tot]=val[i];
        }
        work(1,n,1);
        for(int i=0;i<n;++i) printf("%lld ",h[sth[1]+i]);
    }
}
signed main()
{
    red::main();
    return 0;
}

原文地址:https://www.cnblogs.com/knife-rose/p/12142017.html

时间: 2024-10-09 10:21:21

多项式快速插值的相关文章

luoguP5245 【模板】多项式快速幂

$B(x)=A^k(x)$ $\Rightarrow \ln B(x)=\ln A^k(x)$ $\Rightarrow \ln B(x)=k \ln A(x)$ $\Rightarrow B(x)=\exp(k \ln A(x))$ code: #include <cmath> #include <cstring> #include <algorithm> #include <cstdio> #include <string> #define

@算法 - [email&#160;protected] 多项式的多点求值与快速插值

目录 @0 - 参考资料@ @1 - 多点求值@ @理论推导@ @参考代码@ @例题与应用@ @2 - 快速插值@ @理论推导@ @(不建议参考的)代码@ @例题与应用@(暂无) @0 - 参考资料@ Cyhlnj 的博客 @1 - 多点求值@ @理论推导@ 假设已知多项式 \(A(x)\),使用 FFT 可以将 \(A(w_n^0)\),\(A(w_n^1)\),...,\(A(w_n^{n-1})\) 的值在 \(O(n\log n)\) 的时间内快速求出. 那么问题来了,假如我现在要求解任

多点求值与快速插值

多点求值 给出 $n$ 次多项式 $A(x)$ ,求出 $m$ 个 $x_i$ 对应的 $A(x_i)$ 考虑分治,设 $L(x)=\prod_{i=1}^{\frac{n}{2}}(x-x_i)$ , $R(x)=\prod_{i=\frac{n}{2}+1}^n(x-x_i)$ 对于 $i \in [1,\frac{n}{2}],F(x_i)=(F \mod L)(x_i)$ , 对于 $i \in (\frac{n}{2},n],F(x_i)=(F \mod R)(x_i)$ 就是对于左半

luoguP5219 无聊的水题 I 多项式快速幂

有一个幼儿园容斥:最大次数恰好为 $m=$  最大次数最多为 $m$ - 最大次数最多为 $m-1$. 然后来一个多项式快速幂就好了. code: #include <cmath> #include <cstring> #include <algorithm> #include <cstdio> #include <string> #define ll long long #define ull unsigned long long using

【BZOJ3992】[SDOI2015]序列统计 NTT+多项式快速幂

[BZOJ3992][SDOI2015]序列统计 Description 小C有一个集合S,里面的元素都是小于M的非负整数.他用程序编写了一个数列生成器,可以生成一个长度为N的数列,数列中的每个数都属于集合S. 小C用这个生成器生成了许多这样的数列.但是小C有一个问题需要你的帮助:给定整数x,求所有可以生成出的,且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个.小C认为,两个数列{Ai}和{Bi}不同,当且仅当至少存在一个整数i,满足Ai≠Bi.另外,小C认为这个问题的答案可能

多项式函数插值:多项式形式函数求值的Horner嵌套算法

设代数式序列 $q_1(t), q_2(t), ..., q_{n-1}(t)$ ,由它们生成的多项式形式的表达式(不一定是多项式): $$p(t)=x_1+x_2q_1(t)+...x_nq_1(t)q_2(t)..q_{n-1}(t)=\sum\limits_{i=1}^n(x_i\prod\limits_{j=1}^{i-1}q_j(t))$$ 一般来讲,按照这个形式计算函数在 $t_0$ 点的取值的复杂度为:n-1次 $q_i(t)$ 求值,n-1次浮点数乘法(生成n个不同的乘积),n-

多项式 - 快速沃尔什变换

若\(·\)是一种适用于整数域的二元运算,则两多项式关于此运算的方式定义为 \(C_k = \sum_{i·j=k} A_i * B_j\),即 \(C=A·B\). \(FWT\) 主要解决多项式的常见的三种二元位运算,在三种运算下分别构造出不同的变换方式,个人认为比 \(NTT\) 简单 好背 一些.形式与 \(NTT\) 近似. 没有新东西可说,直接放上洛谷模板题的代码好了: #include <cmath> #include <queue> #include <cst

BZOJ3992 [SDOI2015]序列统计 【生成函数 + 多项式快速幂】

题目 小C有一个集合S,里面的元素都是小于M的非负整数.他用程序编写了一个数列生成器,可以生成一个长度为N的数 列,数列中的每个数都属于集合S.小C用这个生成器生成了许多这样的数列.但是小C有一个问题需要你的帮助: 给定整数x,求所有可以生成出的,且满足数列中所有数的乘积mod M的值等于x的不同的数列的有多少个.小C认为 ,两个数列{Ai}和{Bi}不同,当且仅当至少存在一个整数i,满足Ai≠Bi.另外,小C认为这个问题的答案可能很大 ,因此他只需要你帮助他求出答案mod 1004535809

Luogu5245 【模板】多项式快速幂(多项式exp)

A(x)k=eklnA(x),丝毫不懂为什么指数要对p取模,只是写下exp板子. #include<iostream> #include<cstdio> #include<cmath> #include<cstdlib> #include<cstring> #include<algorithm> using namespace std; #define ll long long #define N 600010 #define P 9