BACKUP
0%

一点简单的多项式常数优化方式

Posted on Edited on
Symbols count in article: 2.1k

简单的看了一下优化多项式常数的东西：source，实现了一下第一部分，这部分比较简单

大体的思路：

  • 两个精度为 \(n\) 的多项式相乘时，可以预先尝试推导最终结果的前 \(n\) 项系数，改变式子形态消除贡献，做长 \(n\) 的循环卷积
  • 消除贡献后去掉多余无贡献部分
  • 将 \(f_0/f_0\) 或 \(f_0/f\) 一类的形态改变成 \(f_0g_0\)，其中 \(g_0\equiv 1/f_0\pmod{x^{n}}\)

以下用无下标的代表模 \(x^{2n}\) 意义下，带下标 \(0\) 的代表模 \(x^n\) 意义下

用 \(\textsf E(n)\) 表示一次长度为 \(n\) 的 FFT 所需时间，例如两个精度为 \(n\) 的多项式相乘需要 \(6\textsf E(n)\)，做长度为 \(2n\) 的循环卷积也需要 \(6\textsf E(n)\)

求逆与除法

即 \(h\equiv f/g\)，直接做需要对 \(g\) 求长 \(2n\) 的逆，已经达到了 \(24\textsf{E}(n)\)，重新审视递推式：

\[ \begin{aligned} h-h_0&\equiv 0&\pmod{x^n} \\ h^2-2hh_0+h_0^2&\equiv 0&\pmod{x^{2n}} \\ fh-2fh_0+gh_0^2&\equiv 0&\pmod{x^{2n}} \\ h&\equiv h_0-\left(\dfrac{g}{f}h_0^2-h_0\right) &\pmod{x^{2n}} \end{aligned} \]

我们希望将 \(f\) 用 \(h_0\) 中的 \(f_0\) 处理掉，注意进行长 \(2n\) 的卷积时，\(h_0\equiv f_0/g_0+\Delta_1 \pmod{x^{2n}}\)，其中 \({\rm ord}(\Delta_1)\geq n\)，因此：

\[ \begin{aligned} \dfrac{g}{f}h_0^2-h_0 \equiv \left(\dfrac{f_0}{f}gh_0 -f_0\right)\dfrac{1}{g_0} +\left(1-\dfrac{g}{f}h_0\right)\Delta_1 &\pmod{x^{2n}} \end{aligned} \]

而 \({\rm ord}(1-(g/f)h_0)\geq n\)，因此 \(\Delta_1\) 不会做任何贡献

\(f_0/f\) 的后 \(n\) 项也会参与贡献，不妨令 \(f_0/f\equiv 1+\Delta_2\)，\(f_1\) 代表 \(f-f_0\)：

\[ \begin{cases} (f_0+f_1)\dfrac{1}{f}\equiv 1 \\ \dfrac{f_0}{f}\equiv 1+\Delta_2 \end{cases} \quad \Rightarrow\quad \Delta_2\equiv-\dfrac{f_1}{f} \]

\[ \dfrac{f_0}{f}gh_0 -f_0\equiv gh_0-f_0-\dfrac{f_1}{f}gh_0 \]

因为 \({\rm ord}(f_1)\geq n\)，\(gh_0/f\) 的精度只需要达到 \(n\) 就可以了，这一部分显然为 \(1\)，最终：

\[ h\equiv h_0-(gh_0-f)\dfrac{1}{g_0} \pmod{x^{2n}} \]

  • 用长 \(n\) 的多项式乘法计算 \(h_0\)：\(6\textsf{E}(n)\)
  • 用长 \(2n\) 的循环卷积计算 \(gh_0-f\)：\(12\textsf{E}(n)\);
  • 长 \(n\) 的求逆计算 \(1/g_0\)：\(12\textsf{E}(n)\)

故模 \(x^n\) 意义下计算 \(h\) 需要 \(15\textsf{E(n)}\)

在倒数的计算过程中，\(h_0\) 和 \(1/g_0 \pmod{x^{n}}\) 是一回事，于是我们不再需要计算 \(h_0\)，同时每次迭代完成后求出的就是 \(1/g_0\pmod{x^{n}}\)，因此从模 \(x^n\) 到模 \(2^n\) 我们仅需要 \(12\textsf{E}(n)\)

但是求逆我们需要从模 \(x^2\) 开始迭代，于是总时间依然是 \(12\textsf{E}(n)\)

exp

求 \(F\) 满足 \(G\circ F\equiv 0\) 我们想让 \(\sum a_n(F-F_0)^n\) 贡献出 \(G\circ F\)，那就将其在 \(F_0\) 处泰勒展开，考虑到 \({\rm ord}(F-F_0)\geq n\)，那么 \((F-F_0)\) 在平方项中就已经不可能向 \(F\) 贡献，即：

\[ \begin{aligned} &G\circ F\equiv G\circ F_0+(G^\prime \circ F_0)(F-F_0)\equiv 0 \\ \Rightarrow\quad&(G^\prime \circ F_0)F_0-G\circ F_0\equiv (G^\prime \circ F_0)F \\ \Rightarrow\quad&F\equiv F_0-\dfrac{G\circ F_0}{G^\prime \circ F_0} \end{aligned} \]

\(g\equiv \exp f\Rightarrow \ln g-f\equiv 0\)，那么 \(G(x)\equiv \ln x-f,G ^\prime\equiv 1/x\)

此时：

\[ g\equiv g_0-g_0\int\left( \dfrac{g_0 ^\prime}{g_0} -f ^\prime\right) \]

直接做的话，我们要求 \(g ^\prime_0/g_0\) 要有 \(2n\) 的精度，于是我们要对 \(g_0\) 求长 \(2n\) 的逆，而这已经达到了 \(24\textsf{E}(n)\)

我们依然期望能将 \(g ^\prime_0/g_0\) 转化成 \(g_0 ^\prime h_0\)，其中 \(h_0\equiv 1/g_0 \pmod{x^n}\)

于是观察 \({\rm ord}\)：\({\rm ord}(g_0h_0 -1)\geq n,{\rm ord}(g_0 ^\prime/g_0-f ^\prime)\geq n-1\)，于是对积分号内部分乘 \(g_0h_0\) 不会有任何影响，于是：

\[ \begin{aligned} \dfrac{g_0 ^\prime}{g_0} -f ^\prime \equiv \dfrac{g_0}{g_0}g_0 ^\prime h_0-g_0h_0 f ^\prime\equiv g_0 ^\prime h_0-f ^\prime-(g_0h_0-1)f_0 ^\prime \end{aligned} \]

其中 \(g_0/g_0\) 的精度达到 \(2n\) 视为 \(1\)，后边 \(f-f_0\) 不会有贡献所以写成 \(f_0\)

  • 用长 \(n\) 的循环卷积计算 \(g_0h_0-1,g_0 ^\prime h_0-f ^\prime\)：\(6\textsf E(n)\)
  • 用长 \(2n\) 的循环卷积计算剩余两次乘法：\(12\textsf{E}(n)\)
  • 用长为 \(n\) 的求逆计算 \(h_0\)：\(12\textsf{E}(n)\)

注意到 \(h_0\) 可以随着 \(\exp\) 的迭代而迭代，且 \(h_0\) 最后一次不需要迭代，因此总时间花费为 \(12+6+6=24\textsf E(n)\)

开根

\[ g\equiv g_0-\dfrac{(g_0^2-f)h_0}{2}\pmod{x^{2n}} \]

  • 用长 \(n\) 的求逆计算 \(h_0\)：\(12\textsf E(n)\)
  • 用长为 \(n\) 的循环卷积计算 \(g_0^2-f\)：\(3\textsf{E}(n)\)
  • 用长为 \(2n\) 的循环卷积计算 \((g_0^2-f)h_0\)：\(6\textsf{E}(n)\)

故花费 \(21\textsf{E}(n)\)，依然注意到 \(h_0\) 也可以随着这个过程迭代，并且最后一轮对 \(h_0\) 的迭代可以省去，求逆过程常数变为原来的 \(1/2\)，故共需要 \(15\textsf E(n)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
namespace Poly {
	template<typename T> void DIF(T *const f,const T *const pr,const int x) {
		for(int len=x,hl=len>>1;hl;len=hl,hl>>=1)
			for(auto s=f,w=pr;s!=f+x;s+=len,w++)
				for(auto l=s,r=s+hl;l!=s+hl;l++,r++) {
					const T tmp=*w**r;
					*r=*l-tmp;
					*l=*l+tmp;
				}
	}
	template<typename T> void DIT(T *const f,const T *const pr,const int x) {
		for(int len=2,hl=1;len<=x;hl=len,len<<=1)
			for(auto s=f,w=pr;s!=f+x;s+=len,w++)
				for(auto l=s,r=s+hl;l!=s+hl;l++,r++) {
					const T ps=*l+*r,ng=*l-*r;
					*l=ps;
					*r=*w*ng;
				}
	}

	template<typename T,int pr> struct poly {
		vector<T> v;
		static T w[(1<<21)+1];

		poly() {}
		poly(T x) {v.push_back(x);}

		static void pre(const int x=1<<21) {
			w[0]=1,w[x]=pr;
			do(i,x>>1,,>>=1) w[i]=w[i<<1]*w[i<<1];
			do(i,x-1) w[i]=w[i&(i-1)]*w[i&-i];
		}

		T operator [] (const int p) const {return p<v.size()?v[p]:T(0);}
		T& operator [] (const int p) {if(p>=len()) v.resize(p+1);return v[p];}

		poly operator + (const poly& rhs) const {return ano()+=rhs;}
		poly operator - (const poly& rhs) const {return ano()-=rhs;}
		template<typename U> poly operator * (const U x) const {return ano()*=x;}
		poly operator / (const poly& rhs) const {return ano()/=rhs;}

		poly& operator - () {
			for(auto &p: v) p=-p;
			return *this;
		}
		poly& operator += (const poly& rhs) {
			if(rhs.len()>len()) v.resize(rhs.len());
			do(i,0,deg()) v[i]+=rhs[i];
			return *this;
		}
		poly& operator -= (const poly& rhs) {
			if(rhs.len()>len()) v.resize(rhs.len());
			do(i,0,deg()) v[i]-=rhs[i];
			return *this;
		}
		template<typename U> poly& operator *= (const U x) {
			for(auto &p: v) p*=x;
			return *this;
		}
		poly& operator *= (const poly& rhs) {
			int n=deg()+rhs.deg(),x=1;
			for(;x<=n;x<<=1);
			return fix_conv(rhs,x);
		}
		poly& operator /= (const int x) {return *this*=T(x).inv();}
		poly& operator /= (const poly& rhs) {
			int n=deg(),x=1;
			for(;x<=n;x<<=1);
			poly sf,sg,sfg;
			sg.cp(rhs,x>>1);
			sfg=sf.cp(*this,x>>1).fix_conv(sg.inv(),x).cut(x>>1);
			sf.fix_conv(rhs,x);
			memset(&sf[0],0,sizeof(T)*(x>>1));
			do(i,x>>1,x-1) sf[i]-=(*this)[i];
			sf.fix_conv(sg,x);
			do(i,0,(x>>1)-1) v[i]=sfg[i];
			do(i,x>>1,x-1) v[i]=-sf[i];
			return *this;
		}

		int deg() const {return v.size()-1;}
		int len() const {return v.size();}
		T* data() {return v.data();}
		poly ano() const {return *this;}
		poly& cp(const poly& rhs,int l) {
			cut(l);
			memcpy(&v[0],&rhs.v[0],sizeof(T)*l);
			return *this;
		}
		poly& cut(const int x) {v.resize(x);return *this;}

		poly& dot(const poly& rhs) {
			if(rhs.len()>len()) v.resize(rhs.len());
			do(i,0,deg()) v[i]*=rhs[i];
			return *this;
		}
		poly& wn(const int x) {
			v.resize(x);
			DIF(data(),w,x);
			return *this;
		}
		poly& revert() {
			DIT(data(),w,len());
			T invx=T(len()).inv();
			do(i,0,deg()) v[i]*=invx;
			reverse(v.begin()+1,v.end());
			return *this;
		}
		poly& fix_conv(poly rhs,const int x) {
			return wn(x).dot(rhs.wn(x)).revert();
		}
		void inv_iter(poly& iv,poly& ib,const int x) const {
			poly ivc=iv;
			ivc.cp(iv,x>>1);
			ib.wn(x).dot(iv.wn(x)).revert();
			memset(&ib[0],0,sizeof(T)*(x>>1));
			ib.wn(x).dot(iv).revert();
			iv.v.swap(ivc.v);
			do(i,x>>1,x-1) iv[i]-=ib[i];
		}
		poly& inv() {
			int n=deg(),x=1;
			for(;x<=n;x<<=1);
			poly g(v[0].inv()),buf;
			do(len,2,<=x,<<=1) inv_iter(g,buf.cp(*this,len),len);
			return *this=g;
		}
		poly& der() {
			do(i,0,deg()-1) v[i]=v[i+1]*(i+1);
			v[deg()]=0;
			return *this;
		}
		poly& intg() {
			do(i,deg(),>=1,--) v[i]=v[i-1]*T(i).inv();
			v[0]=0;
			return *this;
		}
		poly& ln() {return *this=(ano().der()/(*this)).intg();}
		poly& exp() {
			int n=deg(),x=1;
			for(;x<=n;x<<=1);
			poly res(1),gph,gh,giv(1),sf,buf;
			do(len,2,<=x,<<=1) {
				gph=gh=res;
				sf.cp(*this,len).der();
				gph.der().fix_conv(giv,len>>1).cut(len);
				do(i,0,(len>>1)-1) gph[i]-=sf[i]+sf[i+(len>>1)];
				memcpy(&gph[len>>1],&gph[0],sizeof(T)*((len>>1)-1));
				memset(&gph[0],0,sizeof(T)*((len>>1)-1));
				(gh.fix_conv(giv,len>>1))[0]--;
				gh.cut(len);
				memcpy(&gh[len>>1],&gh[0],sizeof(T)*(len>>1));
				memset(&gh[0],0,sizeof(T)*(len>>1));
				gh.fix_conv(sf.cut((len>>1)-1),len);
				memset(&gh[0],0,sizeof(T)*(len>>1));
				(gph-=gh).intg().fix_conv(res,len);
				do(i,len>>1,len-1) res[i]-=gph[i];
				if(len==x) break;
				inv_iter(giv,buf.cp(res,len),len);
			}
			return *this=res;
		}
		poly& sqrt() {
			int n=deg(),x=1;
			for(;x<=n;x<<=1);
			poly res(1),buf,iv(1),ib;
			const T i2=T(2).inv();
			do(len,2,<=x,<<=1) {
				buf=res;
				buf.fix_conv(res,len>>1).cut(len);
				do(i,0,(len>>1)-1) buf[i]-=v[i]+v[i+(len>>1)];
				memcpy(&buf[len>>1],&buf[0],sizeof(T)*(len>>1));
				memset(&buf[0],0,sizeof(T)*(len>>1));
				buf.fix_conv(iv,len);
				do(i,len>>1,len-1) res[i]=-buf[i]*i2;
				if(len==x) break;
				inv_iter(iv,ib.cp(res,len),len);
			}
			return *this=res;
		}
	};

	template<typename T,int pr> T poly<T,pr>::w[(1<<21)+1];

	using gf=poly<mi,31>;
} using namespace Poly;