zoukankan      html  css  js  c++  java
  • 洛谷P4726 【模板】多项式指数函数

    https://www.luogu.org/problemnew/show/P4726

    多项式牛顿迭代http://blog.miskcoo.com/2015/06/polynomial-with-newton-method

    公式:$F(z) equiv F_0(z) - frac{G(F_0(z))}{G'(F_0(z))} pmod {z^n}$

    (以下截图自金策2015论文)

    (研究了很久都没有明白那个泰勒展开的求导为什么是对的...好像是把A(x)当成了”常数项“...?把f看成”变量“???但是g这个函数应该是输入一个函数输出一个函数吧,根本就不是一般的函数,真的能求导吗?但是这种做法的确具有通用性,在多项式开根等也可以用,很迷)

    upd190323:大概感受了一下,以下是口胡:

    首先,假设x是常数。设y=f(x),现在f没有确定,因此y不是常数。g(y)=ln(y)-A(x),其中A(x)是常数

    现在,对于这个确定的x可以用牛顿迭代求出y的零点。$g(y_0)+g'(y_0)(y-y_0)=0$,$y=y_0-frac{g(y_0)}{g'(y_0)}$

    可以发现,对于每一个x,每一次迭代求出来的y是固定的。每一次迭代求出的y与x之间存在一个函数关系。

    因此,如果x不是常数,仍然可以用这样的方法做多项式牛顿迭代。

    不过,就算求导是错的,最后的结论式子也应当是对的?直接两边求ln?(好像也不对)不懂...先鸽着吧)

    式子可以写成$f(x)=f_0(x)+f_0(x)(A(x)-ln(f_0(x))$,可以先算$A(x)-ln(f_0(x))$,根据牛顿迭代那套理论这个后面一半全是0,可以用来简单优化常数

    当然,输入多项式要求常数项是0(因为要求模意义下常数的exp,这个只能对0求),输出多项式常数项就是1

    (式子里面$ln(f_0(x))$求的是长度为len的ln(设$f_0$的长度为len/2,空的部分补0),不是特别懂

    版本1:基于版本3

      1 #prag
      2 ma GCC optimize(2)
      3 #include<cstdio>
      4 #include<algorithm>
      5 #include<cstring>
      6 #include<vector>
      7 #include<cmath>
      8 using namespace std;
      9 #define fi first
     10 #define se second
     11 #define mp make_pair
     12 #define pb push_back
     13 typedef long long ll;
     14 typedef unsigned long long ull;
     15 const int md=998244353;
     16 const int N=262144;
     17 #define delto(a,b) ((a)-=(b),((a)<0)&&((a)+=md))
     18 inline int del(int a,int b)
     19 {
     20     a-=b;
     21     return a<0?a+md:a;
     22 }
     23 int rev[N];
     24 void init(int len)
     25 {
     26     int bit=0,i;
     27     while((1<<(bit+1))<=len)    ++bit;
     28     for(i=0;i<len;++i)
     29         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
     30 }
     31 ull poww(ull a,ull b)
     32 {
     33     ull base=a,ans=1;
     34     for(;b;b>>=1,base=base*base%md)
     35         if(b&1)
     36             ans=ans*base%md;
     37     return ans;
     38 }
     39 void dft(int *a,int len,int idx)//要求len为2的幂
     40 {
     41     int i,j,k,t1,t2;ull wn,wnk;
     42     for(i=0;i<len;++i)
     43         if(i<rev[i])
     44             swap(a[i],a[rev[i]]);
     45     for(i=1;i<len;i<<=1)
     46     {
     47         wn=poww(idx==1?3:332748118,(md-1)/(i<<1));
     48         for(j=0;j<len;j+=(i<<1))
     49         {
     50             wnk=1;
     51             for(k=j;k<j+i;++k,wnk=wnk*wn%md)
     52             {
     53                 t1=a[k];t2=a[k+i]*wnk%md;
     54                 a[k]+=t2;
     55                 (a[k]>=md) && (a[k]-=md);
     56                 a[k+i]=t1-t2;
     57                 (a[k+i]<0) && (a[k+i]+=md);
     58             }
     59         }
     60     }
     61     if(idx==-1)
     62     {
     63         ull ilen=poww(len,md-2);
     64         for(i=0;i<len;++i)
     65             a[i]=a[i]*ilen%md;
     66     }
     67 }
     68 void p_inv(int *f,int *g,int len)//g=f^(-1);f,g数组的长度不小于2len(需要足够长用于临时存放元素) ;要求len是2的幂
     69 {
     70     static int t1[N],t2[N];
     71     g[0]=poww(f[0],md-2);
     72     for(int i=2,j;i<(len<<1);i<<=1)
     73     {
     74         memcpy(t1,f,sizeof(int)*i);
     75         memcpy(t2,g,sizeof(int)*(i>>1));
     76         memset(t2+(i>>1),0,sizeof(int)*(i>>1));
     77         init(i);
     78         dft(t1,i,1);dft(t2,i,1);
     79         for(j=0;j<i;++j)
     80             t1[j]=ull(t1[j])*t2[j]%md;
     81         dft(t1,i,-1);
     82         for(j=0;j<(i>>1);++j)
     83             t1[j]=t1[j+(i>>1)];
     84         memset(t1+(i>>1),0,sizeof(int)*(i>>1));
     85         dft(t1,i,1);
     86         for(j=0;j<i;++j)
     87             t1[j]=ull(t1[j])*t2[j]%md;
     88         dft(t1,i,-1);
     89         for(j=i>>1;j<i;++j)
     90             g[j]=md-t1[j-(i>>1)];
     91     }
     92 }
     93 int inv[300011];
     94 inline void p_de(int *f,int len)//derivative求导;f=f'
     95 {
     96     for(int i=0;i<len-1;++i)
     97         f[i]=ll(i+1)*f[i+1]%md;
     98     f[len-1]=0;
     99 }
    100 inline void p_in(int *f,int len)//integral积分;f=?f
    101 {
    102     for(int i=len-1;i>=1;--i)
    103         f[i]=ll(f[i-1])*inv[i]%md;
    104     f[0]=0;
    105 }
    106 void p_ln(int *f,int len)//要求len为2的幂
    107 {
    108     static int t3[N];
    109     p_inv(f,t3,len);p_de(f,len);
    110     init(len<<1);
    111     dft(f,len<<1,1);dft(t3,len<<1,1);
    112     for(int i=0;i<(len<<1);++i)
    113         f[i]=ull(f[i])*t3[i]%md;
    114     dft(f,len<<1,-1);p_in(f,len);
    115 }
    116 void p_exp(int *f,int *g,int len)//要求len为2的幂,f[0]=0
    117 {
    118     static int t1[N],t2[N],t3[N];
    119     g[0]=1;
    120     for(int i=2,j;i<(len<<1);i<<=1)
    121     {
    122         memcpy(t3,g,sizeof(int)*(i>>1));
    123         memcpy(t2,g,sizeof(int)*(i>>1));
    124         memset(t2+(i>>1),0,sizeof(int)*(i>>1));
    125         memset(g+(i>>1),0,sizeof(int)*(i>>1));
    126         p_ln(g,i);
    127         for(j=0;j<(i>>1);++j)
    128             t1[j]=del(f[j+(i>>1)],g[j+(i>>1)]);
    129         memset(t1+(i>>1),0,sizeof(int)*(i>>1));
    130         init(i);
    131         dft(t1,i,1);dft(t2,i,1);
    132         for(j=0;j<i;++j)
    133             t1[j]=ull(t1[j])*t2[j]%md;
    134         dft(t1,i,-1);
    135         memcpy(g,t3,sizeof(int)*(i>>1));
    136         for(j=(i>>1);j<i;++j)
    137             g[j]=t1[j-(i>>1)];
    138     }
    139 }
    140 int a[N<<1],b[N<<1];
    141 int n,n1;
    142 int main()
    143 {
    144     int i,t;
    145     inv[1]=1;
    146     for(i=2;i<=300000;++i)
    147         inv[i]=ull(md-md/i)*inv[md%i]%md;
    148     scanf("%d",&n);n1=n;
    149     for(i=0;i<n;++i)
    150         scanf("%d",a+i);
    151     for(t=1;t<n;t<<=1);
    152     n=t;
    153     p_exp(a,b,n);
    154     for(i=0;i<n1;++i)
    155         printf("%d ",b[i]);
    156     return 0;
    157 }
    View Code

    版本2:基于此题版本1

      1 #prag
      2 ma GCC optimize(2)
      3 #include<cstdio>
      4 #include<algorithm>
      5 #include<cstring>
      6 #include<vector>
      7 #include<cmath>
      8 using namespace std;
      9 #define fi first
     10 #define se second
     11 #define mp make_pair
     12 #define pb push_back
     13 typedef long long ll;
     14 typedef unsigned long long ull;
     15 const int md=998244353;
     16 const int N=262144;
     17 #define delto(a,b) ((a)-=(b),((a)<0)&&((a)+=md))
     18 inline int del(int a,int b)
     19 {
     20     a-=b;
     21     return a<0?a+md:a;
     22 }
     23 int rev[N];
     24 void init(int len)
     25 {
     26     int bit=0,i;
     27     while((1<<(bit+1))<=len)    ++bit;
     28     for(i=1;i<len;++i)
     29         rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
     30 }
     31 ull poww(ull a,ull b)
     32 {
     33     ull ans=1;
     34     for(;b;b>>=1,a=a*a%md)
     35         if(b&1)
     36             ans=ans*a%md;
     37     return ans;
     38 }
     39 int inv[300011];
     40 void dft(int *a,int len,int idx)//要求len为2的幂
     41 {
     42     int i,j,k,t1,t2;ull wn,wnk;
     43     for(i=0;i<len;++i)
     44         if(i<rev[i])
     45             swap(a[i],a[rev[i]]);
     46     for(i=1;i<len;i<<=1)
     47     {
     48         wn=poww(idx==1?3:332748118,(md-1)/(i<<1));
     49         for(j=0;j<len;j+=(i<<1))
     50         {
     51             wnk=1;
     52             for(k=j;k<j+i;++k,wnk=wnk*wn%md)
     53             {
     54                 t1=a[k];t2=a[k+i]*wnk%md;
     55                 a[k]+=t2;
     56                 (a[k]>=md)&&(a[k]-=md);
     57                 a[k+i]=t1-t2;
     58                 (a[k+i]<0)&&(a[k+i]+=md);
     59             }
     60         }
     61     }
     62     if(idx==-1)
     63     {
     64         ull ilen=inv[len];
     65         for(i=0;i<len;++i)
     66             a[i]=a[i]*ilen%md;
     67     }
     68 }
     69 void p_inv(int *f,int *g,int len)//g=f^(-1);f,g数组的长度不小于2len(需要足够长用于临时存放元素);要求len是2的幂
     70 {
     71     static int t1[N],t2[N];
     72     g[0]=poww(f[0],md-2);
     73     for(int i=2,j;i<=len;i<<=1)
     74     {
     75         memcpy(t1,f,sizeof(int)*i);
     76         memcpy(t2,g,sizeof(int)*(i>>1));
     77         memset(t2+(i>>1),0,sizeof(int)*(i>>1));
     78         init(i);
     79         dft(t1,i,1);dft(t2,i,1);
     80         for(j=0;j<i;++j)
     81             t1[j]=ull(t1[j])*t2[j]%md;
     82         dft(t1,i,-1);
     83         for(j=0;j<(i>>1);++j)
     84             t1[j]=t1[j+(i>>1)];
     85         memset(t1+(i>>1),0,sizeof(int)*(i>>1));
     86         dft(t1,i,1);
     87         for(j=0;j<i;++j)
     88             t1[j]=ull(t1[j])*t2[j]%md;
     89         dft(t1,i,-1);
     90         for(j=i>>1;j<i;++j)
     91             g[j]=md-t1[j-(i>>1)];
     92     }
     93 }
     94 inline void p_de(int *f,int len)//derivative求导;f=f'
     95 {
     96     for(int i=0;i<len-1;++i)
     97         f[i]=ull(i+1)*f[i+1]%md;
     98     f[len-1]=0;
     99 }
    100 inline void p_in(int *f,int len)//integral积分;f=?f
    101 {
    102     for(int i=len-1;i>=1;--i)
    103         f[i]=ull(f[i-1])*inv[i]%md;
    104     f[0]=0;
    105 }
    106 void p_ln(int *f,int len)//要求len为2的幂,f[0]=1
    107 {
    108     static int t3[N];
    109     p_inv(f,t3,len);p_de(f,len);
    110     init(len<<1);
    111     dft(f,len<<1,1);dft(t3,len<<1,1);
    112     for(int i=0;i<(len<<1);++i)
    113         f[i]=ull(f[i])*t3[i]%md;
    114     dft(f,len<<1,-1);p_in(f,len);
    115 }
    116 void p_exp(int *f,int *g,int len)//要求len为2的幂,f[0]=0
    117 {
    118     static int t1[N],t2[N];
    119     g[0]=1;
    120     for(int i=2,j;i<=len;i<<=1)
    121     {
    122         memcpy(t1,g,sizeof(int)*(i>>1));
    123         memset(t1+(i>>1),0,sizeof(int)*(i>>1));
    124         p_ln(t1,i);
    125         for(j=0;j<(i>>1);++j)
    126             t1[j]=del(f[j+(i>>1)],t1[j+(i>>1)]);
    127         memset(t1+(i>>1),0,sizeof(int)*(i>>1));
    128         init(i);
    129         dft(t1,i,1);
    130         memcpy(t2,g,sizeof(int)*(i>>1));
    131         memset(t2+(i>>1),0,sizeof(int)*(i>>1));
    132         dft(t2,i,1);
    133         for(j=0;j<i;++j)
    134             t1[j]=ull(t1[j])*t2[j]%md;
    135         dft(t1,i,-1);
    136         for(j=i>>1;j<i;++j)
    137             g[j]=t1[j-(i>>1)];
    138     }
    139 }
    140 int a[N],b[N];
    141 int n,n1;
    142 int main()
    143 {
    144     int i,t;
    145     inv[1]=1;
    146     for(i=2;i<=300000;++i)
    147         inv[i]=ull(md-md/i)*inv[md%i]%md;
    148     scanf("%d",&n);n1=n;
    149     for(i=0;i<n;++i)
    150         scanf("%d",a+i);
    151     for(t=1;t<n;t<<=1);
    152     n=t;
    153     p_exp(a,b,n);
    154     for(i=0;i<n1;++i)
    155         printf("%d ",b[i]);
    156     return 0;
    157 }
    View Code
  • 相关阅读:
    Java HttpClient使用小结
    【剑指offer】Q18:树的子结构
    poj3041-Asteroids , 二分图的最小顶点覆盖数 = 最大匹配数
    jquery.validate+jquery.form提交的三种方式
    &quot;undefined reference to&quot; 问题解决方法
    [Oracle]
    Effective_java之二:慎用重载函数
    C99规范
    迭代、递归替代循环
    1)Linux程序设计入门--基础知识
  • 原文地址:https://www.cnblogs.com/hehe54321/p/10561943.html
Copyright © 2011-2022 走看看