【学习笔记】多项式求幂/开根

形式化来讲,我们要解决这两个问题:

然而实际上后一个咕了,我曾把这玩意儿当作多项式牛顿迭代的课后题来着……并且更重要的原因没人出这玩意儿的题,所以直接简化成$G^2\equiv F$好了。。。

$1$ 求幂

首先思考式子

怎么化开。那么喜闻乐见的就是一步求对数然后exp回去:

对应的是LuoguP5245:

1
2
3
4
5
6
7
8
9
10
11
signed main(){
int i, y ;
cin >> N, M = 1 ;
Gi = expow(Gp, P - 2) ;
while (M <= N) M <<= 1 ;
cin >> In ; y = strlen(In) ;
for (i = 0 ; i < N ; ++ i) scanf("%lld", &F[i]) ;
for (i = 0 ; i < y ; ++ i) K = (K * 10 + In[i] - '0') % P ;
_Ln(F, L, M) ; for (i = 0 ; i < N ; ++ i) L[i] = 1ll * L[i] * K % P ;
_Exp(L, G, M) ; for (i = 0 ; i < N ; ++ i) printf("%lld ", G[i]) ; return 0 ;
}

但是这种方法只能应用于$a_0=1$的情况下,原因在写exp的时候整理过。

但这东西显然是可以$\Theta(n\log^2n)$暴力ntt做的,我很快乐地水过了这道题的加强版P5273,即不保证$a_0=1$的版本:

似乎大家都是诡异的写法,没有人用$O(n\log^2 n)$去暴力艹这道题。

然而事实上是完全可以卡过去的。我的提交虽然加了-O2#pragma显得十分弟弟,但是其实去掉之后也是很快的。不吹,绝对比大部分的普通$O(n\log^2n)$的NTT跑得快。

那么以下是几个优化的措施:
$\#1$

预处理原根的次幂——卡常利器。

1
2
3
4
5
for (i = 0 ; i < 19  ;++ i){
rr int *rua = gg[i] ; rua[0] = 1 ;
rr int gi = rua[1] = expow(3, 998244352/(1 << (i + 1))) ;
for (j = 2 ; j < (1 << i) ; ++ j) rua[j] = 1ll * rua[j - 1] * gi % P ;
}

然后每次NTT就不需要重新再计算了。

$\#2$

做快速幂的时候注意可以少几次NTT。这点常数优化也是需要的。

1
2
3
4
5
6
7
8
9
10
11
while (K){
NTT(F, M, 1) ;
if (K & 1){
NTT(A, M, 1) ;
for (i = 0 ; i < M ; ++ i) A[i] = 1ll * A[i] * F[i] % P ;
NTT(A, M, -1) ; for (i = N ; i < M ; ++ i) A[i] = 0 ;
}
for (i = 0 ; i < M ; ++ i)
F[i] = 1ll * F[i] * F[i] % P ; NTT(F, M, -1) ;
for (i = N ; i < M ; ++ i) F[i] = 0 ; K >>= 1 ;
}

$\#3$

不用long long.

这其实一个通用的技巧,因为long long一般都比int多好多常数。同时不要强制类型转换而选择1ll*这种形式。实测可以快许多。

不得不说,这个预处理原根确实是一个卡常黑科技qwq

$2$ 开根

普通的开根还是expln

然后就可以胡搞LuoguP5205

1
2
3
4
5
6
void _Invsqr(int *f, int *g, int len){
int Len = 1, l = 0, i, o ;
while (Len < (len << 1)) Len <<= 1, ++ l ; o = expow(2ll, P - 2) ;
for (i = 0 ; i < Len ; ++ i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << l - 1) ;
_Ln(f, E, M) ; for (i = 0 ; i < Len ; ++ i) E[i] = E[i] * o % P ; _Exp(E, g, M) ;
}

但是还是有问题,就是因为需要$\exp$,上述也必须要$F_0=1$。于是当$F_0\neq 1$时,我们可以直接暴力求二次剩余,注意这种方法就需要用牛迭来实现开根

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
namespace solveroot{
#define fr first
#define sc second
#define pint pair<int, int>
#define mp(a, b) make_pair(a, b)
struct complex{ int x, y ; } ;
complex mul(complex a, complex b, int p, int rt){
int nx = (1ll * a.x * b.x % p + 1ll * a.y * b.y % p * rt % p) ;
int ny = (1ll * a.x * b.y % p + 1ll * a.y * b.x % p) ;
return (complex){nx % p, ny % p} ;
}
complex expow(complex a, int b, int p, int z){
complex res = (complex){1, 0} ;
while (b){
if (b & 1)
res = mul(res, a, p, z), b -- ;
a = mul(a, a, p, z), b >>= 1 ;
}
return res ;
}
int expow(int a, int b, int p){
int res = 1 ;
while (b){
if (b & 1)
res = 1ll * a * res % p ;
a = 1ll * a * a % p ; b >>= 1 ;
}
return res ;
}
void seed(){
srand((int)'x'+'q'+'I'+'L'+'o'+'v'+'e'+'u') ;
}
pint solve(int n, int p){
seed() ;
if (!n) return mp(0, 0) ;
int x, y ;
while (1){
x = rand() % p ;
y = (1ll * x * x + (p - n)) % p ;
if (expow(y, (p - 1) >> 1, p) == p - 1) break ;
}
complex now = (complex){x, 1}, ansp, ansn ;
ansp = expow(now, (p + 1) >> 1, p, y),
ansn = (complex){p - ansp.x, ansp.y} ;
if (ansp.x > ansn.x) swap(ansp, ansn) ;
return mp(ansp.x, ansn.x) ;
}
}
。。。。。
LL Ig[MAXN], pf[MAXN] ;
using namespace solveroot ;
void _sqr(LL *f, LL *g, int len){
if (len <= 1){
g[0] = solve(f[0], P).fr ;
return ;
}
rr int i, l = 0, Len = 1 ;
_sqr(f, g, (len + 1) >> 1) ;
while (Len < (len << 1)) Len <<= 1, ++ l ;
for (i = len ; i < Len ; ++ i) pf[i] = Ig[i] = 0 ;
for (i = 0 ; i < len ; ++ i) Ig[i] = 0, pf[i] = 2 * g[i] % P ;
_Inv(pf, Ig, len) ; for (i = 0 ; i < Len ; ++ i) /* */ t[i] = 0 ;
for (i = 1 ; i < Len ; ++ i) R[i] = (R[i >> 1] >> 1) | ((i & 1) << ( l - 1 )) ;
NTT(g, Len, 1) ; for (i = 0 ; i < Len ; ++ i) g[i] = g[i] * g[i] % P ;
NTT(g, Len, -1) ; for (i = 0 ; i < len ; ++ i) g[i] = (g[i] + f[i]) % P ;
NTT(g, Len, 1) ; NTT(Ig, Len, 1) ; for (i = 0 ; i < Len ; ++ i) g[i] = g[i] * Ig[i] % P ;
NTT(g, Len, -1) ; for (i = len ; i < Len ; ++ i) g[i] = 0 ;
}

然后为了做这道题还去学了Cipolla,然后还顺便发现luogu上少一道二次剩余的板子就顺手出了道题2333