• 洛谷P4726 【模板】多项式指数函数


    传送门

    orz这题太珂怕了……似乎都找不到几个板子参(chao)考(xi)……前置芝士又特别多……而且我写的时候牛顿迭代那里NTT数组长度写错了调了半天……然后各种地方多项式没清零又调了半天……可能是因为平时都抄板子的缘故没注意这问题……

    前置芝士:多项式对数函数(这里),泰勒展开(可以看看这里第一个回答,非常……生动形象),牛顿迭代法(这里

    泰勒展开

    简单来说的话,就是要求一个函数$f(x)$某一点上的值,我们可以构造一个函数$g(x)$,那么对于$x$点的值,有$$f(x)approx g(x)=g(0)+frac{f^1(0)}{1!}x+frac{f^2(0)}{2!}x^2+……+frac{f^n(0)}{n!}x^n$$

    其中$f^n(0)$表示对原函数图像上$0$这个点进行$n$阶求导

    牛顿迭代

    牛顿迭代可以用来求一个函数的零点,多项式牛顿迭代自然是用来求多项式的零点的,即对于一个函数$G(x)$,求满足条件$G(F(z)) equiv 0 pmod {z^n}$的多项式$F(z)$

    miskcoo大佬是这么说的

    扯远了

    然后现在就是要计算$$F(x)=e^{A(x)}$$

    那么变形一下得$$ln F(x)-A(x)=0$$

    我们设$G(F(x))=ln F(x)-A(x)$,那么就是要求这一个函数的零点。那么我们把$F(x)$看做变量,$A(x)$看做常数(我也不知道为什么能这样),对这个进行求导,得$G'(F(x))=frac{1}{F(x)}$

    那么代入上面牛顿迭代的公式得$$F(x)equiv F_0(x)-frac{G(F_0(x))}{G'(F_0(x))}pmod{x^n}$$
    $$F(x)equiv F_0(x)(1-ln F_0(x)+A(x))pmod{x^n}$$

    然后因为$A(0)=0$,所以$F(x)$的常数项为1

    然后左转把各种板子复制过来就好了

      1 //minamoto
      2 #include<iostream>
      3 #include<cstdio>
      4 #include<algorithm>
      5 #define swap(x,y) (x^=y,y^=x,x^=y)
      6 #define mul(x,y) (1ll*x*y%P)
      7 #define add(x,y) (x+y>=P?x+y-P:x+y)
      8 #define dec(x,y) (x-y<0?x-y+P:x-y)
      9 using namespace std;
     10 #define getc() (p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++)
     11 char buf[1<<21],*p1=buf,*p2=buf;
     12 inline int read(){
     13     #define num ch-'0'
     14     char ch;bool flag=0;int res;
     15     while(!isdigit(ch=getc()))
     16     (ch=='-')&&(flag=true);
     17     for(res=num;isdigit(ch=getc());res=res*10+num);
     18     (flag)&&(res=-res);
     19     #undef num
     20     return res;
     21 }
     22 char sr[1<<21],z[20];int K=-1,Z;
     23 inline void Ot(){fwrite(sr,1,K+1,stdout),K=-1;}
     24 inline void print(int x){
     25     if(K>1<<20)Ot();if(x<0)sr[++K]=45,x=-x;
     26     while(z[++Z]=x%10+48,x/=10);
     27     while(sr[++K]=z[Z],--Z);sr[++K]=' ';
     28 }
     29 const int N=500005,P=998244353;
     30 inline int ksm(int a,int b){
     31     int res=1;
     32     while(b){
     33         if(b&1) res=mul(res,a);
     34         a=mul(a,a),b>>=1;
     35     }
     36     return res;
     37 }
     38 int n,r[N],A[N],B[N],C[N],D[N],F[N],G[N],O[N],f[N],g[N],inv[N];
     39 void init(int limit){
     40     for(int i=0,l=limit<<1;i<=l;++i) inv[i]=ksm(i,P-2);
     41 }
     42 void NTT(int *A,int type,int len){
     43     int limit=1,l=0;
     44     while(limit<len) limit<<=1,++l;
     45     for(int i=0;i<limit;++i)
     46     r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
     47     for(int i=0;i<limit;++i)
     48     if(i<r[i]) swap(A[i],A[r[i]]);
     49     for(int mid=1;mid<limit;mid<<=1){
     50         int R=mid<<1,Wn=ksm(3,(P-1)/R);O[0]=1;
     51         for(int j=1;j<mid;++j) O[j]=mul(O[j-1],Wn);
     52         for(int j=0;j<limit;j+=R){
     53             for(int k=0;k<mid;++k){
     54                 int x=A[j+k],y=mul(O[k],A[j+k+mid]);
     55                 A[j+k]=add(x,y),A[j+k+mid]=dec(x,y);
     56             }
     57         }
     58     }
     59     if(type==-1){
     60         reverse(A+1,A+limit);
     61         for(int i=0,invl=inv[limit];i<limit;++i)
     62         A[i]=mul(A[i],invl);
     63     }
     64 }
     65 void Inv(int *a,int *b,int len){
     66     if(len==1) return (void)(b[0]=inv[a[0]]);
     67     Inv(a,b,len>>1);
     68     for(int i=0;i<len;++i) C[i]=a[i],D[i]=b[i];
     69     NTT(C,1,len<<1),NTT(D,1,len<<1);
     70     for(int i=0,l=(len<<1);i<l;++i) C[i]=mul(mul(C[i],D[i]),D[i]);
     71     NTT(C,-1,len<<1);
     72     for(int i=0;i<len;++i) b[i]=dec(add(b[i],b[i]),C[i]);
     73     for(int i=0,l=(len<<1);i<l;++i) C[i]=D[i]=0;
     74 }
     75 void Direv(int *A,int *B,int len){
     76     for(int i=1;i<len;++i) B[i-1]=mul(A[i],i);B[len-1]=0; 
     77 }
     78 void Inter(int *A,int *B,int len){
     79     for(int i=1;i<len;++i) B[i]=mul(A[i-1],inv[i]);B[0]=0; 
     80 }
     81 void Ln(int *a,int *b,int len){
     82     Direv(a,A,len),Inv(a,B,len);
     83     NTT(A,1,len<<1),NTT(B,1,len<<1);
     84     for(int i=0,l=len<<1;i<l;++i) A[i]=mul(A[i],B[i]);
     85     NTT(A,-1,len<<1),Inter(A,b,len<<1);
     86     for(int i=0,l=len<<1;i<l;++i) A[i]=B[i]=0;
     87 }
     88 void Exp(int *a,int *b,int len){
     89     if(len==1) return (void)(b[0]=1);
     90     Exp(a,b,len>>1),Ln(b,F,len);
     91     F[0]=dec(a[0]+1,F[0]);
     92     for(int i=1;i<len;++i) F[i]=dec(a[i],F[i]);
     93     NTT(F,1,len<<1),NTT(b,1,len<<1);
     94     for(int i=0,l=len<<1;i<l;++i) b[i]=mul(b[i],F[i]);
     95     NTT(b,-1,len<<1);
     96     for(int i=len,l=(len<<1);i<l;++i) b[i]=F[i]=0;
     97 }
     98 int main(){
     99 //    freopen("testdata.in","r",stdin);
    100     n=read();
    101     for(int i=0;i<n;++i) f[i]=read();
    102     int len;for(len=1;len<=n;len<<=1);init(len);
    103     Exp(f,g,len);
    104     for(int i=0;i<n;++i) print(g[i]);
    105     Ot();
    106     return 0;
    107 }
  • 相关阅读:
    sql: 生日三个月内有效
    asp 日期操作
    Csharp:操作存儲過程輸出參數,和返回值
    sql: 生日赠品中的相关算法
    asp and javascript: sql server export data to csv and to xls
    Csharp:user WebControl Read Adobe PDF Files In Your Web Browser
    Csharp: Send Email
    数据结构(七)排序---冒泡排序
    数据结构(七)排序---基本概念和分类
    数据结构(六)查找---散列表(哈希表)查找
  • 原文地址:https://www.cnblogs.com/bztMinamoto/p/9749384.html
Copyright © 2020-2023  润新知