首页 > 代码库 > 多项式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
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。