首页 > 代码库 > 多项式exp

多项式exp

写了一发多项式exp,picks太神辣

  1 #include <bits/stdc++.h>
  2 
  3 using namespace std;
  4 
  5 typedef long long ll;
  6 const int N = 262144,mod = 998244353,G = 3;
  7 
  8 void read(int&x){
  9     x = 0;char ch = getchar();
 10     while (ch < 0 || 9 < ch) ch = getchar();
 11     while (0 <= ch&& ch <= 9) x = x * 10 + ch - 0,ch = getchar(); 
 12 }
 13 int power(int x,int y,int z){
 14     int ans = 1;
 15     while (y){
 16         if (y&1) ans = (ll)ans*x%z;
 17         x = (ll)x*x%z;
 18         y >>= 1;
 19     }
 20     return ans;
 21 }
 22 int add(int x,int y){
 23     x += y;
 24     while (x >= mod) x -= mod;
 25     return x;
 26 }
 27 
 28 
 29 int a[N+10],b[N+10],c[N+10],d[N+10],p[N];
 30 int pg[N+10],pg_inv[N+10],inv[N+10];
 31 int n,m;
 32 void FFT(int*,int,int);
 33 void poly_exp(int*,int*,int*,int*,int);
 34 void poly_inv(int*,int*,int*,int);
 35 void poly_ln(int*,int*,int*,int);
 36 void poly_der(int*,int);
 37 void poly_mot(int*,int); 
 38 int main(){
 39     read(n);
 40     for (int i = 0;i < n;i++) read(a[n-i-1]);
 41     for (m = 1;m < n;m<<=1);
 42     for(int i = 1;i <= (m<<1);i <<= 1){
 43         pg[i] = power(G,(mod-1)/i,mod);
 44         pg_inv[i] = power(pg[i],mod-2,mod);
 45     }
 46     inv[1] = 1;
 47     for (int i = 2;i <= N;i++)inv[i] = mod-(ll)mod/i*inv[mod%i]%mod;
 48     /* 
 49     poly_inv(a,b,c,m); 
 50 //    for (int i = 0;i < m<<1;i++)cout << a[i]<<‘ ‘;cout<<endl;
 51 //    for (int i = 0;i < m<<1;i++)cout << b[i]<<‘ ‘;cout<<endl;
 52 //    cout <<(ll)a[0]*b[0]%mod<<endl;
 53     FFT(a,1,m<<1);
 54 //    for (int i = 0;i < m<<1;i++)cout <<a[i]<<‘*‘;cout<<endl;
 55     FFT(b,1,m<<1);
 56     for(int i = 0;i <m<<1;i++) a[i] =(ll)a[i]*b[i]%mod;
 57     FFT(a,-1,m<<1);
 58 //    for (int i = 0;i < m;i++)cout << a[i]<<‘ ‘;cout<<endl;
 59     bool flag = 1;
 60     for(int i = 1;i < m;i++)if (a[i] != 0) flag = 0;
 61     if(a[0] != 1) flag = 0;
 62     puts(flag ? "YES" : "NO");
 63     */ 
 64 //    /*
 65     poly_exp(a,b,c,d,m);
 66 //    for (int i = 0;i < m<<1;i++) cout << a[i]<<‘ ‘;cout << endl;
 67 //    for (int i = 0;i < m<<1;i++)cout << b[i] <<‘ ‘;cout<<endl;
 68     /*
 69 //    b[0] = 1;b[1] = 1;b[2] = inv[2];b[3] = inv[6];
 70     poly_ln(b,c,d,m);bool flag = 1;
 71 //    for (int i = 0;i < m;i++) cout << c[i]<<‘ ‘;cout<<endl;
 72     for (int i = 0;i <m;i++)if(c[i] != a[i]) flag = 0;
 73     puts(flag ? "YES" : "NO");
 74     */
 75 //    */
 76     return 0;
 77 }
 78 void FFT(int *a,int d,int deg){
 79     for (int i = 1;i < deg;i <<= 1)
 80         for (int j = 0;j < i;j++)
 81             p[i+j] = p[j]+deg/i/2;
 82     for (int i = 0;i <deg;i++) if (p[i] > i) swap(a[i],a[p[i]]);
 83     for (int n = 2;n <= deg;n <<= 1){
 84         int wn = (d == 1 ? pg[n] : pg_inv[n]);
 85         for (int i = 0;i <deg;i += n){
 86             int w = 1;
 87             for (int j = i;j < i + n/2;j++){
 88                 int x = a[j],y = (ll)a[j+n/2]*w%mod;
 89                 a[j] = add(x,y);
 90                 a[j+n/2] = add(x,mod-y);
 91                 w = (ll)w*wn%mod;
 92             }
 93         }
 94     }
 95     if (d == -1){
 96         for (int i = 0;i < deg;i++)
 97             a[i] = (ll)a[i]*inv[deg]%mod;
 98     }
 99 }
100 void poly_der(int*a,int deg){
101     for (int i = 0;i < deg;i++)
102         a[i] = (ll)a[i+1]*(i+1)%mod;
103 }
104 void poly_mot(int*a,int deg){
105     for (int i = deg;i > 0;i--)
106         a[i] = (ll)a[i-1]*inv[i]%mod;
107     a[0] = 0;
108 }
109 void poly_inv(int*a,int*b,int*c,int deg){
110     if (deg == 1){
111         if (a[0] <= N) b[0] = inv[a[0]];else 
112         b[0] = power(a[0],mod-2,mod);
113         return;
114     }
115     poly_inv(a,b,c,deg/2);
116     int deg2 = deg<<1;
117     //B = B(2-BA)
118     for (int i = 0;i < deg;i++) c[i] = a[i];
119     for (int i = deg;i < deg2;i++) c[i] = 0;
120     FFT(c,1,deg2);
121     FFT(b,1,deg2);
122     for (int i = 0;i < deg2;i++) b[i] = (ll)b[i]*add(2,mod-(ll)b[i]*c[i]%mod)%mod;
123     FFT(b,-1,deg2);
124     for (int i = deg;i < deg2;i++) b[i] = 0;
125 }
126 void poly_ln(int*a,int*b,int*c,int deg){
127     /*
128         b = lna
129         b‘ = a‘/a
130     */
131     int deg2 = deg<<1;
132     for (int i = 0;i < deg2;i++)c[i] = 0;
133     poly_inv(a,c,b,deg);
134     for(int i = 0;i <= deg;i++) b[i] = a[i];
135     for(int i = deg+1;i < deg2;i++) b[i] = 0;
136     poly_der(b,deg);
137     FFT(b,1,deg2);
138     FFT(c,1,deg2);
139     for (int i = 0;i < deg2;i++) b[i] = (ll)b[i]*c[i]%mod;
140     FFT(b,-1,deg2);
141     for (int i = deg;i < deg2;i++) b[i] = 0;
142     poly_mot(b,deg);
143     for(int i = deg;i < deg2;i++)b[i] = 0;
144 }
145 void poly_exp(int *a,int *b,int *c,int*d,int deg){
146     //B = B(1-lnB+A)
147     if(deg == 1){
148         //if a[0] = 0
149         b[0] = 1;
150         return;
151     }
152     int deg2=deg<<1;    
153     poly_exp(a,b,c,d,deg/2);
154     for (int i = deg;i < deg2;i++) b[i] = 0;
155     poly_ln(b,c,d,deg);
156 
157     for (int i = 0;i < deg;i++) c[i] = add(a[i]+(i==0),mod-c[i]);
158 //    cout <<"1-lnB+A ";for (int i = 0;i <deg2;i++) cout << c[i]<<‘ ‘;cout<<endl;    
159     for (int i = deg;i < deg2;i++) c[i] = 0;
160     FFT(c,1,deg2);
161     FFT(b,1,deg2);
162     for(int i = 0;i < deg2;i++) b[i] = (ll)b[i]*c[i]%mod;
163     FFT(b,-1,deg2);
164     for (int i = deg;i < deg2;i++) b[i] = 0;
165 //    for (int i = 0;i < deg2;i++)cout << b[i]<<‘ ‘;cout<<"deg:"<<deg<<"\n";
166 }

因为只会求e^0的值,所以只能解决常数项为0的多项式的exp。。。。

多项式exp其实求的是e^A的泰勒展开的前若干项,

用牛顿迭代法,求出解的递推式技术分享

最后一行在多项式的意义下在这个情况中证明了这个方法至少是平方收敛(多项式下的牛顿迭代应该都是可以类似证精度的)(数意义下的不会证orz)

然后。。我试图用原始的初等的推多项式求逆的那一套方法推这个式子,成功GG。。。。

然后写的时候发现一些问题,就是这个求原函数啊,求模x^n意义下的原函数需要模x^(n+1)意义下的函数,这题因为需要求原函数的函数都是没取模的所以不影响,不然还有注意一下细节。。。

多项式exp