参考文献：Jörg Arndt，Matter s Computational，DOI 10.1007/978-3-642-14764-7
整数运算
Karatsuba algorithm
对于两个
n
n
n比特整数
A
=
x
a
1
+
a
0
B
=
x
b
1
+
b
0
\begin{aligned} A = xa_1+a_0\\ B = xb_1+b_0 \end{aligned}
A=xa1+a0B=xb1+b0
乘积为
A
B
=
(
x
2
+
x
)
a
1
b
1
+
(
x
+
1
)
a
0
b
0
+
x
(
a
1
−
a
0
)
(
b
0
−
b
1
)
AB = (x^2+x)a_1b_1 + (x+1)a_0b_0 + x(a_1-a_0)(b_0-b_1)
AB=(x2+x)a1b1+(x+1)a0b0+x(a1−a0)(b0−b1)
将
n
n
n比特数字相乘，转化为3个
n
/
2
n/2
n/2比特数字相乘。因此，它的复杂度为
O
(
n
log
3
)
≈
O
(
n
1.585
)
O(n^{\log3}) \approx O(n^{1.585})
O(nlog3)≈O(n1.585)
这个算法是 2-way splitting 的，即数字被分为 a 1 , a 0 a_1,\,a_0 a1,a0两部分。
Toom-Cook algorithm，3-way splitting，复杂度为 O ( n log 3 5 ) ≈ O ( n 1.465 ) O(n^{\log_3 5}) \approx O(n^{1.465}) O(nlog35)≈O(n1.465)
Bodrato-Zanoni algorithm，4-way splitting，复杂度为 O ( n log 4 7 ) ≈ O ( n 1.403 ) O(n^{\log_4 7}) \approx O(n^{1.403}) O(nlog47)≈O(n1.403)
FFT、NTT
数字也是多项式。数字 a a a的 R R R进制表示： a n − 1 a n − 2 ⋯ a 0 = ∑ i = 0 n − 1 a i R i a_{n-1}a_{n-2} \cdots a_0 = \sum_{i=0}^{n-1} a_i R^i an−1an−2⋯a0=∑i=0n−1aiRi
我们利用NTT (无精度损失)，在频域上计算阵列乘，等价于在时域上计算卷积。
复杂度为： 2 O ( n log n ) + O ( n ) = O ( n log n ) 2\,O(n \log n)+O(n) = O(n \log n) 2O(nlogn)+O(n)=O(nlogn)
在使用NTT时，对于 Z p [ x ] / ( x 2 n − 1 ) Z_p[x]/(x^{2n}-1) Zp[x]/(x2n−1)， 2 n ∣ p − 1 2n \mid p-1 2n∣p−1，要注意保证 p > n ( R − 1 ) 2 p > n(R-1)^2 p>n(R−1)2，使得计算结果正确。如果要连续多次乘法， p p p需要取非常大的值才行。
如果使用FFT，在计算过程中是复数运算，不可避免的损失精度。但是计算机中double类型可以表示很大范围的数字，所以没有 p p p过大无法表示的烦恼。
Power
若
e
=
e
s
⋯
e
1
e
0
e=e_s \cdots e_1 e_0
e=es⋯e1e0，那么
a
e
=
(
a
)
e
0
(
a
2
)
e
1
⋯
(
a
2
s
)
e
s
a^e = (a)^{e_0} (a^2)^{e_1} \cdots (a^{2^s})^{e_s}
ae=(a)e0(a2)e1⋯(a2s)es
Right-to-Left算法：从
e
0
e_0
e0向
e
s
e_s
es扫描，并不断计算
(
a
2
j
)
2
(a^{2^j})^2
(a2j)2；扫到了
e
j
=
1
e_j=1
ej=1，需要计算
r
a
2
j
ra^{2^j}
ra2j
Left-to-Right算法：从 e s e_s es向 e 0 e_0 e0扫描，并不断计算 r = r 2 r=r^2 r=r2；扫到了 e j = 1 e_j=1 ej=1，需要计算 r a ra ra
如果 a a a的规模较小，那么后者会比前者快数倍！
这里的整数乘积和平方可以使用FFT加速；由于Left-to-Right中的 a a a不变，因此其FFT只需要做一次 (FFT-caching)，之后的乘法都在FFT域进行。
浮点数运算
Inverse
对于高精度浮点数除法，代价过于昂贵。可以使用浮点数乘法，迭代逼近它的逆。
给定
d
d
d，设置初始值
x
0
≈
1
d
x_0 \approx \dfrac{1}{d}
x0≈d1 (两三个有效位的逼近即可)，迭代
x
k
+
1
=
x
k
(
1
+
(
1
−
d
x
k
)
)
x_{k+1} = x_k(1+(1-dx_k))
xk+1=xk(1+(1−dxk))
注意，不要化简！计算机中有精度损失，会出现相抵消现象。
令 ∣ e ∣ < 1 |e|<1 ∣e∣<1是相对错误大小。如果 x k = 1 + e d x_k = \dfrac{1+e}{d} xk=d1+e，那么 x k + 1 = 1 − e 2 d x_{k+1} = \dfrac{1-e^2}{d} xk+1=d1−e2
随着迭代进行，相对错误指数级减小。
Inverse square root
给定
d
d
d，设置初始值
x
0
≈
1
d
x_0 \approx \dfrac{1}{\sqrt d}
x0≈d1 ，迭代
x
k
+
1
=
x
k
(
1
+
1
−
d
x
k
2
2
)
x_{k+1} = x_k(1 + \frac{1-dx_k^2}{2})
xk+1=xk(1+21−dxk2)
令
e
e
e是相对错误大小。如果
x
k
=
1
+
e
d
x_k = \dfrac{1+e}{\sqrt d}
xk=d1+e，那么
x
k
+
1
=
1
−
3
e
2
/
2
−
e
3
/
2
d
x_{k+1} = \dfrac{1-3e^2/2 - e^3/2}{\sqrt d}
xk+1=d1−3e2/2−e3/2
Inverse Cube root
给定
d
d
d，设置初始值
x
0
≈
1
d
1
/
3
x_0 \approx \dfrac{1}{d^{1/3}}
x0≈d1/31 ，迭代
x
k
+
1
=
x
k
(
1
+
1
−
d
2
x
k
3
3
)
x_{k+1} = x_k(1 + \frac{1-d^2 x_k^3}{3})
xk+1=xk(1+31−d2xk3)
令
e
e
e是相对错误大小。如果
x
k
=
1
+
e
d
1
/
3
x_k = \dfrac{1+e}{d^{1/3}}
xk=d1/31+e，那么
x
k
+
1
=
1
−
2
e
2
/
2
−
4
e
3
/
3
−
e
4
/
3
d
1
/
3
x_{k+1} = \dfrac{1-2e^2/2 - 4e^3/3 - e^4/3}{d^{1/3}}
xk+1=d1/31−2e2/2−4e3/3−e4/3
扩展到Group
一、定义函数：
ϕ
(
x
)
:
=
x
(
1
+
(
1
−
d
x
)
)
\phi(x) := x(1+(1-dx))
ϕ(x):=x(1+(1−dx))
假设已知
x
0
d
=
1
m
o
d
p
x_0d = 1 \mod p
x0d=1modp，那么
x
1
=
ϕ
(
x
0
)
=
ϕ
(
d
−
1
(
1
+
k
p
)
)
=
d
−
1
(
1
−
k
2
p
2
)
=
d
−
1
m
o
d
p
2
x_1 = \phi(x_0) = \phi(d^{-1}(1+kp)) = d^{-1}(1-k^2p^2) = d^{-1} \mod p^2
x1=ϕ(x0)=ϕ(d−1(1+kp))=d−1(1−k2p2)=d−1modp2
进一步的，
x
k
=
ϕ
(
x
k
−
1
)
=
d
−
1
m
o
d
p
2
k
x_k = \phi(x_{k-1}) = d^{-1} \mod p^{2^k}
xk=ϕ(xk−1)=d−1modp2k
二、定义函数：
ϕ
(
x
)
:
=
x
(
1
+
2
−
1
(
1
−
d
x
2
)
)
\phi(x) := x(1 + 2^{-1}(1-dx^2))
ϕ(x):=x(1+2−1(1−dx2))
假设已知
x
0
2
d
=
1
m
o
d
p
x_0^2d = 1 \mod p
x02d=1modp，那么
x
1
=
ϕ
(
x
0
)
=
ϕ
(
d
−
1
/
2
(
1
+
k
p
)
)
=
d
−
1
/
2
(
1
−
k
2
p
2
)
=
d
−
1
/
2
m
o
d
p
2
x_1 = \phi(x_0) = \phi(d^{-1/2}(1+kp)) = d^{-1/2}(1-k^2p^2) = d^{-1/2} \mod p^2
x1=ϕ(x0)=ϕ(d−1/2(1+kp))=d−1/2(1−k2p2)=d−1/2modp2
进一步的，
x
k
=
ϕ
(
x
k
−
1
)
=
d
−
1
/
2
m
o
d
p
2
k
x_k = \phi(x_{k-1}) = d^{-1/2} \mod p^{2^k}
xk=ϕ(xk−1)=d−1/2modp2k
实质：级数展开
一、令
y
=
1
−
d
x
y=1-dx
y=1−dx，那么
1
d
=
x
⋅
1
1
−
y
\dfrac{1}{d} = x \cdot \dfrac{1}{1-y}
d1=x⋅1−y1
级数展开，记
ϕ
k
(
x
)
:
=
x
(
1
+
y
+
⋯
+
y
k
−
1
)
\phi_k(x) := x(1+y+\cdots+y^{k-1})
ϕk(x):=x(1+y+⋯+yk−1)
那么
ϕ
k
(
1
+
e
d
)
=
1
−
(
−
e
)
k
d
\phi_k(\dfrac{1+e}{d}) = \dfrac{1 - (-e)^k}{d}
ϕk(d1+e)=d1−(−e)k
易知
ϕ
2
(
x
)
=
x
(
1
+
y
)
=
x
(
1
+
(
1
−
d
x
)
)
\phi_2(x)=x(1+y) = x(1+(1-dx))
ϕ2(x)=x(1+y)=x(1+(1−dx))
二、令
y
=
1
−
d
x
2
y=1-dx^2
y=1−dx2，那么
1
d
=
x
⋅
1
1
−
y
\dfrac{1}{\sqrt d} = x \cdot \dfrac{1}{\sqrt{1-y}}
d1=x⋅1−y1
级数展开，记
ϕ
k
+
1
(
x
)
:
=
x
(
1
+
y
2
+
⋯
+
(
k
2
k
)
y
k
4
k
)
\phi_{k+1}(x) := x(1 + \dfrac{y}{2} + \cdots + \dfrac{(_k^{2k})y^{k}}{4^k})
ϕk+1(x):=x(1+2y+⋯+4k(k2k)yk)
易知
ϕ
2
(
x
)
=
x
(
1
+
y
2
)
=
x
(
1
+
1
−
d
x
2
2
)
\phi_2(x) = x(1+\dfrac{y}{2}) = x(1+\dfrac{1-dx^2}{2})
ϕ2(x)=x(1+2y)=x(1+21−dx2)
r-th root
一、关于
d
1
/
r
d^{1/r}
d1/r的2阶迭代公式
ϕ
2
(
x
)
=
x
+
d
−
x
r
r
x
r
−
1
=
(
r
−
1
)
x
+
d
x
r
−
1
r
\phi_2(x) = x + \dfrac{d-x^r}{rx^{r-1}} = \frac{(r-1)x+\dfrac{d}{x^{r-1}}}{r}
ϕ2(x)=x+rxr−1d−xr=r(r−1)x+xr−1d
若
e
→
0
±
e \rightarrow 0^\pm
e→0±，那么
ϕ
2
(
d
1
/
r
⋅
1
+
e
1
−
e
)
=
d
1
/
r
⋅
(
1
+
e
1
−
e
+
(
1
−
e
)
r
−
(
1
+
e
)
r
r
⋅
(
1
+
e
)
r
−
1
(
1
−
e
)
)
≈
d
1
/
r
⋅
(
1
+
e
1
−
e
+
(
1
−
r
e
)
−
(
1
+
r
e
)
r
⋅
(
1
+
r
e
−
e
)
(
1
−
e
)
)
≈
d
1
/
r
⋅
1
+
(
r
−
3
)
e
+
r
e
2
1
+
(
r
−
3
)
e
+
(
r
−
2
)
e
2
\begin{aligned} \phi_2(d^{1/r} \cdot \dfrac{1+e}{1-e}) &= d^{1/r} \cdot (\dfrac{1+e}{1-e} + \dfrac{(1-e)^r - (1+e)^r} {r \cdot (1+e)^{r-1} (1-e)}) \\ &\approx d^{1/r} \cdot (\dfrac{1+e}{1-e} + \dfrac{(1-re) - (1+re)} {r \cdot (1+re-e) (1-e)}) \\ &\approx d^{1/r} \cdot \dfrac{1+(r-3)e+re^2}{1+(r-3)e+(r-2)e^2}\\ \end{aligned}
ϕ2(d1/r⋅1−e1+e)=d1/r⋅(1−e1+e+r⋅(1+e)r−1(1−e)(1−e)r−(1+e)r)≈d1/r⋅(1−e1+e+r⋅(1+re−e)(1−e)(1−re)−(1+re))≈d1/r⋅1+(r−3)e+(r−2)e21+(r−3)e+re2
二、迭代算法
-
设置 x 0 = d , E 0 = d r − 1 x_0=d,\, E_0=d^{r-1} x0=d,E0=dr−1
-
迭代，直到 x k x_k xk足够逼近 x ∞ = d 1 / r x_\infty = d^{1/r} x∞=d1/r：
P k : = 1 + 1 − E k r → 1 x k + 1 : = x k ⋅ P k E k + 1 : = E K ⋅ P k r → 1 \begin{aligned} P_k &:= 1+\dfrac{1-E_k}{r} &\rightarrow 1\\ x_{k+1} &:= x_k \cdot P_k\\ E_{k+1} &:= E_K \cdot P_k^r &\rightarrow 1\\ \end{aligned} Pkxk+1Ek+1:=1+r1−Ek:=xk⋅Pk:=EK⋅Pkr→1→1 -
最后，得到 x k ≈ d 1 / r x_k \approx d^{1/r} xk≈d1/r
还有类似这种迭代算法的更高阶 (Higher order) 迭代算法，略。
多项式
迭代器的收敛速率
对于多项式 f ( x ) f(x) f(x)，若 r r r是一个根 f ( r ) = 0 f(r)=0 f(r)=0
关于 r r r的**(单点) 迭代器**表示为： x k + 1 = ϕ ( x k ) , x ∞ = r x_{k+1}=\phi(x_k),\, x_\infty = r xk+1=ϕ(xk),x∞=r
它需要满足：fixed point， ϕ ( r ) = r \phi(r)=r ϕ(r)=r；attracting， ∣ ϕ ′ ( r ) ∣ < 1 |\phi'(r)| < 1 ∣ϕ′(r)∣<1
另外，多点迭代器表示为： x k + 1 = ϕ ( x k , x k − 1 , ⋯ , x k − j ) x_{k+1}=\phi(x_k,x_{k-1},\cdots,x_{k-j}) xk+1=ϕ(xk,xk−1,⋯,xk−j)
例如，割线法：
x
k
+
1
=
ϕ
(
x
k
,
x
k
−
1
)
=
x
k
−
f
(
x
k
)
⋅
x
k
−
x
k
−
1
f
(
x
k
)
−
f
(
x
k
−
1
)
x_{k+1}=\phi(x_k,x_{k-1}) = x_k - f(x_k)\cdot \dfrac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}
xk+1=ϕ(xk,xk−1)=xk−f(xk)⋅f(xk)−f(xk−1)xk−xk−1
令
x
=
r
(
1
+
e
)
,
∣
e
∣
≪
1
x=r(1+e),\,|e| \ll 1
x=r(1+e),∣e∣≪1，如果
ϕ
(
x
)
=
r
(
1
+
α
e
n
+
O
(
e
n
+
1
)
)
\phi(x)=r(1+\alpha e^n+O(e^{n+1}))
ϕ(x)=r(1+αen+O(en+1))，那么
ϕ
\phi
ϕ是
n
o
r
d
e
r
n\,\,order
norder迭代器。若
n
=
1
n=1
n=1，则是线性的 (linear)；否则，叫做超线性的 (super-linear)，表现更好。
对于
n
n
n阶迭代器
ϕ
\phi
ϕ，如果
f
(
r
)
=
0
f(r)=0
f(r)=0，那么
r
r
r是super-attracting fixed point：
ϕ
′
(
r
)
=
0
,
ϕ
′
′
(
r
)
=
0
,
⋯
,
ϕ
(
n
−
1
)
(
r
)
=
0
\phi'(r)=0,\,\phi''(r)=0,\,\cdots,\,\phi^{(n-1)}(r)=0
ϕ′(r)=0,ϕ′′(r)=0,⋯,ϕ(n−1)(r)=0
因此，
ϕ
(
x
)
+
f
(
x
)
n
⋅
h
(
x
)
\phi(x) + f(x)^n \cdot h(x)
ϕ(x)+f(x)n⋅h(x)与
ϕ
(
x
)
\phi(x)
ϕ(x)有相同的收敛阶数，
h
(
x
)
h(x)
h(x)是定义在
r
r
r邻域上的任意函数。
Sample root
如果
f
(
x
)
f(x)
f(x)只有单根，令
ϕ
(
x
)
:
=
x
−
p
(
x
)
f
(
x
)
,
p
(
x
)
:
=
f
′
(
x
)
−
1
m
o
d
f
(
x
)
\phi(x) := x-p(x)f(x),\, p(x) := f'(x)^{-1} \mod f(x)
ϕ(x):=x−p(x)f(x),p(x):=f′(x)−1modf(x)
其中
p
(
x
)
p(x)
p(x)是
f
′
(
x
)
f'(x)
f′(x)的模逆，因此
deg
p
(
x
)
<
deg
f
(
x
)
\deg p(x) < \deg f(x)
degp(x)<degf(x)，于是
deg
ϕ
(
x
)
≤
2
deg
f
(
x
)
−
1
\deg \phi(x) \le 2\deg f(x) -1
degϕ(x)≤2degf(x)−1
那么， ϕ ( x ) \phi(x) ϕ(x)是 f ( x ) f(x) f(x)的根的二阶迭代器， x k = ϕ ( x k ) , f ( x ∞ ) = 0 x_k=\phi(x_k),\, f(x_\infty)=0 xk=ϕ(xk),f(x∞)=0
如果
p
f
′
+
q
f
≡
1
pf'+qf \equiv 1
pf′+qf≡1，那么更高阶迭代器为：
p
1
:
=
p
ϕ
1
:
=
x
−
p
1
f
p
r
:
=
p
p
r
−
1
′
−
q
(
r
−
1
)
p
r
−
1
ϕ
r
:
=
ϕ
r
−
1
+
(
−
1
)
r
p
r
f
r
r
!
\begin{aligned} p_1 &:= p\\ \phi_1 &:= x-p_1f\\ p_r &:= pp'_{r-1} - q(r-1)p_{r-1}\\ \phi_r &:= \phi_{r-1} + \dfrac{(-1)^rp_rf^r}{r!} \end{aligned}
p1ϕ1prϕr:=p:=x−p1f:=ppr−1′−q(r−1)pr−1:=ϕr−1+r!(−1)rprfr
其中
ϕ
r
\phi_r
ϕr是
f
(
x
)
f(x)
f(x)的根的
r
−
1
r-1
r−1阶迭代器。
Multiple root
对于重根，迭代器 ϕ ( x ) = x − f f ′ \phi(x) = x-\dfrac{f}{f'} ϕ(x)=x−f′f不收敛。
例如，对于 f ( x ) = ( x 2 − d ) m f(x)=(x^2-d)^m f(x)=(x2−d)m，有 ϕ ( x ) = x − x 2 − d 2 m x \phi(x) = x - \dfrac{x^2-d}{2mx} ϕ(x)=x−2mxx2−d，它只有当 m > 1 m>1 m>1时才会收敛。
m
m
m重根的二阶迭代器为
ϕ
2
(
x
)
:
=
x
−
m
⋅
f
f
′
\phi_2(x) := x - m \cdot \dfrac{f}{f'}
ϕ2(x):=x−m⋅f′f
如果不知道根
r
r
r的重数，用
F
=
f
/
f
′
F=f / f'
F=f/f′替换
f
f
f即可。
其实，
F
F
F与
f
f
f有相同的根集合，但
F
F
F的所有根都是单根：
f
o
r
f
(
r
)
=
0
,
∀
r
l
e
t
f
(
x
)
=
(
x
−
r
)
m
h
(
x
)
h
(
r
)
≠
0
t
h
e
n
F
(
x
)
=
f
(
x
)
/
f
′
(
x
)
=
(
x
−
r
)
h
(
x
)
m
h
(
x
)
+
(
x
−
r
)
h
′
(
x
)
\begin{aligned} for\\ & f(r) &=& 0,\,\, \forall r\\ let\\ & f(x) &=& (x-r)^mh(x)\\ & h(r) &\not =& 0\\ then\\ & F(x) &=& f(x)/f'(x)\\ & &=& (x-r)\dfrac{h(x)}{mh(x)+(x-r)h'(x)} \end{aligned}
forletthenf(r)f(x)h(r)F(x)=====0,∀r(x−r)mh(x)0f(x)/f′(x)(x−r)mh(x)+(x−r)h′(x)h(x)
获得高阶迭代器
令
ϕ
2
(
x
)
=
x
−
f
(
x
)
/
f
′
(
x
)
\phi_2(x) = x - f(x)/f'(x)
ϕ2(x)=x−f(x)/f′(x)是2阶迭代器，那么计算
ϕ
r
(
x
)
=
ϕ
r
−
1
(
x
)
−
f
(
ϕ
r
−
1
(
x
)
)
f
′
(
x
)
\phi_r(x) = \phi_{r-1}(x) - \dfrac{f(\phi_{r-1}(x))}{f'(x)}
ϕr(x)=ϕr−1(x)−f′(x)f(ϕr−1(x))
将得到
r
r
r阶迭代器
ϕ
r
(
x
)
\phi_r(x)
ϕr(x)
扫一扫
494
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
