Processing math: 100%
Yijan's Blog HomeArchives
多项式全家桶

暂时只更了 多项式求逆 ，多项式 Ln ，多项式 exp。。

多项式求逆

众所周知，一个多项式求逆后是一个无穷次项的玩意。但是我们一般只需要它的前 n 项。

我们需要算出：

A×B≡1(modxn)

这东西怎么做？背板 我们可以考虑倍增来求解。

首先如果 A 仅有一项，那么 B[0] 就是 A[0] 的逆元。

否则，假设我们已经有了

A×B′≡1(modx⌈n2⌉)

同时我们必然有真正的 B 在模 x⌊n2⌋ 的时候也是 1 ，也就是说有

A×B≡1(modx⌈n2⌉)

那么考虑两个式子相减

A×(B−B′)≡0(modx⌈n2⌉)

由于 A 显然不同余于 0 ，所以有

B−B′≡0(modx⌈n2⌉)

所以我们可以平方一下，把模的位置翻倍，并且由于是向上取整，我们知道前 n 项都会是 0。

(B−B′)2≡0(modxn)B2−2BB′+B′2≡0(modxn)

然后我们可以给式子的一遍乘上一个 A ，仍然前 n 项是 0 ：

A(B2−2BB′+B′2)≡0(modxn)

由于 AB 显然是 1 ，所以我们拆开有

B−2B′+AB′2≡0(modxn)

所以移项后有

B≡2B′−AB′2(modxn)B≡B′(2−AB′)(modxn)

这个东西复杂度用主定理分析出来是 O(nlogn) 而不是 O(nlog2n) 哦。

同时也正好复习一下 NTT 的板子。等会顺便也复习下 分治 NTT 吧 QwQ

求逆跑得好快啊

cpp
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
#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "cmath"
#include "vector"
#include "map"
#include "set"
#include "queue"
using namespace std;
#define MAXN 300006
//#define int long long
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; ++i)
#define per(i, a, b) for (int i = (a), i##end = (b); i >= i##end; --i)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define vi vector<int>
#define all(x) (x).begin() , (x).end()
#define mem( a ) memset( a , 0 , sizeof a )
#define P 998244353
typedef long long ll;
int n , m;
int A[MAXN] , B[MAXN];

int Pow( int x , int a ) {
    int ans = 1 , cur = x % P;
    while( a ) {
        if( a & 1 ) ans = 1ll * ans * cur % P;
        cur = 1ll * cur * cur % P , a >>= 1;
    }
    return ans;
}

int Wn[2][MAXN] , rev[MAXN];
void getwn( int len ) {
    for( int mid = 1 ; mid < len ; mid <<= 1 ) {
        int w0 = Pow( 3 , ( P - 1 ) / ( mid << 1 ) ) , w1 = Pow( 3 , P - 1 - ( P - 1 ) / ( mid << 1 ) );
        Wn[0][mid] = Wn[1][mid] = 1;
        for( int i = 1 ; i < mid ; ++ i )
            Wn[0][mid + i] = 1ll * Wn[0][mid + i - 1] * w0 % P,
            Wn[1][mid + i] = 1ll * Wn[1][mid + i - 1] * w1 % P;
    }
}
void getr( int len ) {
    int t = __builtin_ctz( len ) - 1;
    for( int i = 1 ; i < len ; ++ i ) rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << t );
}
void NTT( int* A , int len , int typ ) {
    rep( i , 0 , len - 1 ) if( i < rev[i] ) swap( A[i] , A[rev[i]] );
    for( int mid = 1 ; mid < len ; mid <<= 1 )
        for( int i = 0 ; i < len ; i += ( mid << 1 ) )
            for( int k = 0 ; k < mid ; ++ k ) {
                int t1 = A[i + k] , t2 = 1ll * A[i + k + mid] * Wn[typ][mid + k] % P;
                A[i + k] = (t1 + t2) % P , A[i + k + mid] = (t1 + P - t2) % P;
            }
    if( typ == 1 ) for( int inv = Pow( len , P - 2 ) , i = 0 ; i < len ; ++ i ) A[i] = 1ll * A[i] * inv % P;
}

int tmpa[MAXN];
inline void Inv( int* A , int* B , int n ) { // B = inv(A) mod x^n
    if( n == 1 ) return (void) ( B[0] = Pow( A[0] , P - 2 ) );
    Inv( A , B , ( n + 1 ) >> 1 );
    int len = 1 , l = 0;
    while( len <= n * 2 ) len <<= 1 , ++ l;
    getr( len ) , getwn( len );
    rep( i , 0 , n - 1 ) tmpa[i] = A[i];
    rep( i , n , len - 1 ) tmpa[i] = 0;
    NTT( tmpa , len , 0 ) , NTT( B , len , 0 );
    rep( i , 0 , len - 1 ) B[i] = 1ll * B[i] * ( 2 - 1ll * tmpa[i] * B[i] % P + P ) % P;
    NTT( B , len , 1 );
    rep( i , n , len - 1 ) B[i] = 0;
}

void solve() {
    cin >> n;
    rep( i , 0 , n - 1 ) scanf("%d",A + i);
    Inv( A , B , n );
    rep( i , 0 , n - 1 ) printf("%d ",B[i]);
}

signed main() {
//    freopen("input","r",stdin);
//    freopen("fuckout","w",stdout);
//    int T;cin >> T;while( T-- ) solve();
    solve();
}

多项式 Ln

我们考虑

B=lnA

两边导一下：

B′=1AA′

所以我们可以对 A 求逆乘上 A′ 然后再做个积分就完事了。

cpp
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
#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "cmath"
#include "vector"
#include "map"
#include "set"
#include "queue"
using namespace std;
#define MAXN 300006
//#define int long long
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; ++i)
#define per(i, a, b) for (int i = (a), i##end = (b); i >= i##end; --i)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define vi vector<int>
#define all(x) (x).begin() , (x).end()
#define mem( a ) memset( a , 0 , sizeof a )
#define P 998244353
typedef long long ll;
int n , m;
int A[MAXN] , B[MAXN];

int Pow( int x , int a ) {
    int ans = 1 , cur = x % P;
    while( a ) {
        if( a & 1 ) ans = 1ll * ans * cur % P;
        cur = 1ll * cur * cur % P , a >>= 1;
    }
    return ans;
}

int Wn[2][MAXN] , rev[MAXN];
void getwn( int len ) {
    for( int mid = 1 ; mid < len ; mid <<= 1 ) {
        int w0 = Pow( 3 , ( P - 1 ) / ( mid << 1 ) ) , w1 = Pow( 3 , P - 1 - ( P - 1 ) / ( mid << 1 ) );
        Wn[0][mid] = Wn[1][mid] = 1;
        for( int i = 1 ; i < mid ; ++ i )
            Wn[0][mid + i] = 1ll * Wn[0][mid + i - 1] * w0 % P,
            Wn[1][mid + i] = 1ll * Wn[1][mid + i - 1] * w1 % P;
    }
}
void getr( int len ) {
    int t = __builtin_ctz( len ) - 1;
    for( int i = 1 ; i < len ; ++ i ) rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << t );
}
void NTT( int* A , int len , int typ ) {
    rep( i , 0 , len - 1 ) if( i < rev[i] ) swap( A[i] , A[rev[i]] );
    for( int mid = 1 ; mid < len ; mid <<= 1 )
        for( int i = 0 ; i < len ; i += ( mid << 1 ) )
            for( int k = 0 ; k < mid ; ++ k ) {
                int t1 = A[i + k] , t2 = 1ll * A[i + k + mid] * Wn[typ][mid + k] % P;
                A[i + k] = (t1 + t2) % P , A[i + k + mid] = (t1 + P - t2) % P;
            }
    if( typ == 1 ) for( int inv = Pow( len , P - 2 ) , i = 0 ; i < len ; ++ i ) A[i] = 1ll * A[i] * inv % P;
}

int tmpa[MAXN];
inline void Inv( int* A , int* B , int n ) { // B = inv(A) mod x^n
    if( n == 1 ) return (void) ( B[0] = Pow( A[0] , P - 2 ) );
    Inv( A , B , ( n + 1 ) >> 1 );
    int len = 1 , l = 0;
    while( len <= n * 2 ) len <<= 1 , ++ l;
    getr( len ) , getwn( len );
    rep( i , 0 , n - 1 ) tmpa[i] = A[i];
    rep( i , n , len - 1 ) tmpa[i] = 0;
    NTT( tmpa , len , 0 ) , NTT( B , len , 0 );
    rep( i , 0 , len - 1 ) B[i] = 1ll * B[i] * ( 2 - 1ll * tmpa[i] * B[i] % P + P ) % P;
    NTT( B , len , 1 );
    rep( i , n , len - 1 ) B[i] = 0;
}

int ta[MAXN] , invA[MAXN];
inline void Ln( int* A , int* B , int n ) {
    mem( ta ) , mem( invA );
    Inv( A , invA , n );
    rep( i , 0 , n - 1 ) ta[i] = 1ll * A[i + 1] * ( i + 1 ) % P;
    int len = 1 , l = 0;
    while( len <= n + n ) len <<= 1 , ++ l;
    getr( len ) , getwn( len );
    NTT( ta , len , 0 ) , NTT( invA , len , 0 );
    rep( i , 0 , len - 1 ) ta[i] = 1ll * ta[i] * invA[i] % P;
    NTT( ta , len , 1 );
    rep( i , 1 , n - 1 ) B[i] = 1ll * ta[i - 1] * Pow( i , P - 2 ) % P;
}

void solve() {
    cin >> n;
    rep( i , 0 , n - 1 ) scanf("%d",A + i);
    Ln( A , B , n );
    rep( i , 0 , n - 1 ) printf("%d ",B[i]);
}

signed main() {
//    freopen("input","r",stdin);
//    freopen("fuckout","w",stdout);
//    int T;cin >> T;while( T-- ) solve();
    solve();
}

多项式 exp

已知 B(x)≡expA(x)(modxn) ，给定 A(x) 求出 B(x)。

首先考虑构造 G(x) ，满足 G(B(x))=0 ，那么

G(B(x))=lnB(x)−A(x)

我们知道

G′(B(x))=1B(x)

那么带入泰勒展开那里的公式

Bt+1(x)=Bt(x)−G(Bt(x))G′(Bt(x))(modx2t+1)=Bt(x)−lnBt(x)−A(x)1Bt(x)(modx2t+1)=Bt(x)−Bt(x)(lnBt(x)−A(x))(modx2t+1)=Bt(x)[1−lnBt(x)+A(x)](modx2t+1)

实现上，我们可以先递归下去算出 B′(x)≡A(x)(mod⌈n2⌉) ，然后通过上面那个式子把它变成 B(x)(modxn) 。

复杂度也可以分析得 O(nlogn) 但是常数很大。（虽然据说论文哥可以让它跑得比我写的 NTT 还快）

cpp
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
#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "cmath"
#include "vector"
#include "map"
#include "set"
#include "queue"
using namespace std;
#define MAXN 300006
//#define int long long
#define rep(i, a, b) for (int i = (a), i##end = (b); i <= i##end; ++i)
#define per(i, a, b) for (int i = (a), i##end = (b); i >= i##end; --i)
#define pii pair<int,int>
#define fi first
#define se second
#define mp make_pair
#define pb push_back
#define eb emplace_back
#define vi vector<int>
#define all(x) (x).begin() , (x).end()
#define mem( a ) memset( a , 0 , sizeof a )
#define P 998244353
typedef long long ll;
int n , m;
int A[MAXN] , B[MAXN];

int Pow( int x , int a ) {
    int ans = 1 , cur = x % P;
    while( a ) {
        if( a & 1 ) ans = 1ll * ans * cur % P;
        cur = 1ll * cur * cur % P , a >>= 1;
    }
    return ans;
}

int Wn[2][MAXN] , rev[MAXN];
void getwn( int len ) {
    for( int mid = 1 ; mid < len ; mid <<= 1 ) {
        int w0 = Pow( 3 , ( P - 1 ) / ( mid << 1 ) ) , w1 = Pow( 3 , P - 1 - ( P - 1 ) / ( mid << 1 ) );
        Wn[0][mid] = Wn[1][mid] = 1;
        for( int i = 1 ; i < mid ; ++ i )
            Wn[0][mid + i] = 1ll * Wn[0][mid + i - 1] * w0 % P,
                    Wn[1][mid + i] = 1ll * Wn[1][mid + i - 1] * w1 % P;
    }
}
void getr( int len ) {
    int t = __builtin_ctz( len ) - 1;
    for( int i = 1 ; i < len ; ++ i ) rev[i] = ( rev[i >> 1] >> 1 ) | ( ( i & 1 ) << t );
}
void NTT( int* A , int len , int typ ) {
    rep( i , 0 , len - 1 ) if( i < rev[i] ) swap( A[i] , A[rev[i]] );
    for( int mid = 1 ; mid < len ; mid <<= 1 )
        for( int i = 0 ; i < len ; i += ( mid << 1 ) )
            for( int k = 0 ; k < mid ; ++ k ) {
                int t1 = A[i + k] , t2 = 1ll * A[i + k + mid] * Wn[typ][mid + k] % P;
                A[i + k] = (t1 + t2) % P , A[i + k + mid] = (t1 + P - t2) % P;
            }
    if( typ == 1 ) for( int inv = Pow( len , P - 2 ) , i = 0 ; i < len ; ++ i ) A[i] = 1ll * A[i] * inv % P;
}

int tmpa[MAXN];
inline void Inv( int* A , int* B , int n ) { // B = inv(A) mod x^n
    if( n == 1 ) return (void) ( B[0] = Pow( A[0] , P - 2 ) );
    Inv( A , B , ( n + 1 ) >> 1 );
    int len = 1 , l = 0;
    while( len <= n * 2 ) len <<= 1 , ++ l;
    getr( len ) , getwn( len );
    rep( i , 0 , n - 1 ) tmpa[i] = A[i];
    rep( i , n , len - 1 ) tmpa[i] = 0;
    NTT( tmpa , len , 0 ) , NTT( B , len , 0 );
    rep( i , 0 , len - 1 ) B[i] = 1ll * B[i] * ( 2 - 1ll * tmpa[i] * B[i] % P + P ) % P;
    NTT( B , len , 1 );
    rep( i , n , len - 1 ) B[i] = 0;
}

int ta[MAXN] , invA[MAXN];
inline void Ln( int* A , int* B , int n ) {
    mem( ta ) , mem( invA );
    Inv( A , invA , n );
    rep( i , 0 , n - 1 ) ta[i] = 1ll * A[i + 1] * ( i + 1 ) % P;
    int len = 1 , l = 0;
    while( len <= n + n ) len <<= 1 , ++ l;
    getr( len ) , getwn( len );
    NTT( ta , len , 0 ) , NTT( invA , len , 0 );
    rep( i , 0 , len - 1 ) ta[i] = 1ll * ta[i] * invA[i] % P;
    NTT( ta , len , 1 );
    rep( i , 1 , n - 1 ) B[i] = 1ll * ta[i - 1] * Pow( i , P - 2 ) % P;
}

int tln[MAXN];
inline void Exp( int* A , int* B , int n ) { // B equals to exp(A)
    if( n == 1 ) return void( B[0] = 1 );
    Exp( A , B , ( n + 1 ) >> 1 );
    int len = 1 , l = 0;
    while( len <= 2 * n ) len <<= 1 , ++ l;
    for( int i = 0 ; i < len ; ++ i ) tln[i] = 0;
    Ln( B , tln , n );
    rep( i , 0 , n - 1 ) tln[i] = ( A[i] - tln[i] + P ) % P;
    ++ tln[0];
    getr( len ) , getwn( len );
    NTT( B , len , 0 ) , NTT( tln , len , 0 );
    rep( i , 0 , len - 1 ) B[i] = 1ll * B[i] * tln[i] % P;
    NTT( B , len , 1 );
    rep( i , n , len - 1 ) B[i] = 0;
}

void solve() {
    cin >> n;
    rep( i , 0 , n - 1 ) scanf("%d",A + i);
    Exp( A , B , n );
    rep( i , 0 , n - 1 ) printf("%d ",B[i]);
}

signed main() {
//    freopen("input","r",stdin);
//    freopen("fuckout","w",stdout);
//    int T;cin >> T;while( T-- ) solve();
    solve();
}
文章作者: yijan
文章链接: https://yijan.co/duoxinagshi/
版权声明: 本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 Yijan's Blog！
NTT总结多项式求逆