洛谷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 }
版本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 }
