多項式各種操作

前些天學了一堆多項式的算法，今天總結一下。

##乘法

樸素的NTT就不說了，主要說一下三模數NTT。

三模數NTT，顧名思義，就是選取三個适合做NTT的模數，然後把它們用CRT合并起來得到的答案再去對我們要求的模數取模。

為了友善，這三個模數分别是998244353, 1004535809和469762049。它們的原根都是3，且它們減去1的值都有超過20個2的因子。

如果你覺得能用int128的話就請忽略下面所有的表述

但是，在用CRT合并的時候，我們遇到一個麻煩：這三個模數的乘積太大了。我們令三個模數分别為p1, p2, p3，對三個模數取模得到的答案分别為a1, a2, a3，不做取模的原本的答案為Ans，則有：

A n s ≡ a 1   ( m o d   p 1 ) Ans \equiv a1~(mod~p1) Ans≡a1 (mod p1)

A n s ≡ a 2   ( m o d   p 2 ) Ans \equiv a2~(mod~p2) Ans≡a2 (mod p2)

A n s ≡ a 3   ( m o d   p 3 ) Ans \equiv a3~(mod~p3) Ans≡a3 (mod p3)

先用CRT合并前兩個模數的答案，得到

A n s ≡ A   ( m o d   M ) Ans \equiv A~(mod~M) Ans≡A (mod M)

A n s ≡ a 3   ( m o d   p 3 ) Ans \equiv a3~(mod~p3) Ans≡a3 (mod p3)

設

A n s = k M + A Ans = kM + A Ans=kM+A

那麼

k M + A ≡ a 3   ( m o d   p 3 ) kM + A \equiv a3~(mod~p3) kM+A≡a3 (mod p3)

即

k ≡ ( a 3 − A ) M   ( m o d   p 3 ) k \equiv \frac{(a3 - A)}{M}~(mod~p3) k≡M(a3−A)​ (mod p3)

此時k, M, A我們都已經得到，就直接對我們要求的模數取模就好了，即

A n s   %   p = ( ( k   %   p ) ∗ ( M   %   p ) + A )   %   p Ans~\%~p = ((k~\%~p) * (M~\%~p) + A)~\%~p Ans % p=((k % p)∗(M % p)+A) % p

##Code

inline void exntt(int *a, int *b)
{
    int len = 1, bit = 0;
    while (len < n + m) len <<= 1, bit++;
    for (int i = 0; i < len; i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << bit - 1;
    for (int i = 0; i < 3; i++) {
        for (int j = 0; j < len; j++) c[j] = a[j], d[j] = b[j];
        ntt(c, len, 1, p[i]);
        ntt(d, len, 1, p[i]);
        for (int j = 0; j < len; j++) ans[i][j] = 1ll * c[j] * d[j] % p[i];
        ntt(ans[i], len, -1, p[i]);
    }
    for (int i = 0; i < n + m - 1; i++) {
        ll A = (mul(1ll * ans[0][i] * p[1] % M, ksm(p[1] % p[0], p[0] - 2, p[0]), M) +
                mul(1ll * ans[1][i] * p[0] % M, ksm(p[0] % p[1], p[1] - 2, p[1]), M)) % M;
        ll k = ((ans[2][i] - A) % p[2] + p[2]) % p[2] * ksm(M % p[2], p[2] - 2, p[2]) % p[2];
        printf("%lld ", ((k % mod) * (M % mod) % mod + A % mod) % mod);
    }
}
           

##逆元

多項式求逆是個好東西，基本上所有除了乘法外的東西都要用到求逆。

我們要求

A B ≡ 1   m o d   x n AB \equiv 1~mod~x^n AB≡1 mod xn

假如我們已經求得

A B ′ ≡ 1   m o d   x n 2 AB' \equiv 1~mod~x^{\frac{n}{2}} AB′≡1 mod x2n​

則有

B ≡ B ′   m o d   x n 2 B \equiv B'~mod~x^{\frac{n}{2}} B≡B′ mod x2n​

把 B ′ B' B′移到左邊然後平方得到

B 2 − 2 B B ′ + B ′ 2 ≡ 0   m o d   x n B^2 - 2BB' + B'^2 \equiv 0~mod~x^n B2−2BB′+B′2≡0 mod xn

B 2 B^2 B2不好處理，是以同乘以一個A，即

B − 2 B ′ + A B ′ 2 ≡ 0   m o d   x n B - 2B' + AB'^2 \equiv 0~mod~x^n B−2B′+AB′2≡0 mod xn

我們就得到了

B ≡ 2 B ′ − A B ′ 2   m o d   x n B \equiv 2B' - AB'^2~mod~x^n B≡2B′−AB′2 mod xn

于是倍增求即可。

##Code

inline void poly_inv(int *a, int *b, int n)
{
    if (n == 1) {
        b[0] = ksm(a[0], mod - 2);
        return;
    }
    poly_inv(a, b, n + 1 >> 1);
    int len = 1;
    while (len < n << 1) len <<= 1;
    for (int i = 0; i < n; i++) tmp[i] = a[i];
    for (int i = n; i < len; i++) tmp[i] = 0;
    ntt(tmp, len, 1);
    ntt(b, len, 1);
    for (int i = 0; i < len; i++)
        b[i] = 1ll * b[i] * (2 - 1ll * b[i] * tmp[i] % mod + mod) % mod;
    ntt(b, len, -1);
    for (int i = n; i < len; i++) b[i] = 0;
}
           

##開方

多項式開方的思路和求逆一樣，都是倍增的利用上一次的答案求下一次的。

簡單推一下式子吧：

B 2 ≡ A   m o d   x n B^2 \equiv A~mod~x^n B2≡A mod xn

B ′ 2 ≡ A   m o d   x n 2 B'^2 \equiv A~mod~x^{\frac{n}{2}} B′2≡A mod x2n​

B 2 − B ′ 2 ≡ 0   m o d   x n 2 B^2 - B'^2 \equiv 0~mod~x^{\frac{n}{2}} B2−B′2≡0 mod x2n​

這裡出現了兩組解，我們隻保留

B − B ′ ≡ 0   m o d   x n 2 B - B' \equiv 0~mod~x^{\frac{n}{2}} B−B′≡0 mod x2n​

還是套路的平方

B 2 − 2 B B ′ + B ′ 2 ≡ 0   m o d   x n B^2 - 2BB' + B'^2 \equiv 0~mod~x^n B2−2BB′+B′2≡0 mod xn

其實 B 2 B^2 B2就是 A A A

A − 2 B B ′ + B ′ 2 ≡ 0   m o d   x n A - 2BB' + B'^2 \equiv 0~mod~x^n A−2BB′+B′2≡0 mod xn

那麼最後 B B B就求出來了

B ≡ A + B ′ 2 2 B ′   m o d   x n B \equiv \frac{A + B'^2}{2B'}~mod~x^n B≡2B′A+B′2​ mod xn

##Code

void poly_sqrt(int *a, int *b, int n)
{
    if (n == 1) {
        b[0] = 1; //一般情況下a[0] = 1
        return;
    }
    poly_sqrt(a, b, n + 1 >> 1);
    int len = 1;
    while (len < n << 1) len <<= 1;
    memset(c, 0, sizeof c);
    poly_inv(b, c, n);
    for (int i = 0; i < n; i++) tmp[i] = a[i];
    for (int i = n; i < len; i++) tmp[i] = 0;
    ntt(tmp, len, 1);
    ntt(b, len, 1);
    ntt(c, len, 1);
    for (int i = 0; i < len; i++) b[i] = 1ll * (tmp[i] + 1ll * b[i] * b[i] % mod) * c[i] % mod * inv2 % mod;
    ntt(b, len, -1);
    for (int i = n; i < len; i++) b[i] = 0;
}
           

##除法和取模

已經有了求逆，那麼多項式除法和取模又是什麼？和逆元有什麼差別嗎？

舉個例子：

由國小數學得到

157 12 = 13 , 157 ≡ 1   m o d   12 \frac{157}{12} = 13, 157 \equiv 1~mod~12 12157​=13,157≡1 mod 12

然而， 12 ∗ 13 12*13 12∗13并不等于 157 157 157，是以 157 157 157乘以 12 12 12的逆元也并不等于 13 13 13，這就是除法和求逆的差別。

那除法怎麼求呢？

把 157 , 12 , 13 , 1 157, 12, 13, 1 157,12,13,1都看成多項式，我們發現是取模的結果 1 1 1産生了上述影響。

如果我們把被除數多項式（ 157 157 157）、除數多項式（ 12 12 12）和商多項式（ 13 13 13）都翻轉一下，就可以得到

21 ∗ 31 = 651 , 751 − 651 = 100 21*31=651, 751 - 651 = 100 21∗31=651,751−651=100

我們發現，此時取模的結果 1 1 1也被翻轉成了100，而100的位數較高，不會産生導緻除法和求逆結果不同的影響。

是以，我們隻要用 751 751 751去乘上 21 21 21的逆元就可以得到 31 31 31，然後把 31 31 31翻轉回來就得到了 13 13 13。

什麼？你說取模？

你都求出了 13 13 13，你還不會用 157 − 12 ∗ 13 157-12*13 157−12∗13來求 1 1 1嗎？

##Code

inline void poly_div(int *a, int *b, int *d, int *r, int n, int m)
{
    for (int i = 0; i < n; i++) a1[i] = a[n - i - 1];
    for (int i = 0; i < m; i++) b1[i] = b[m - i - 1];
    int l = n - m + 1;
    poly_inv(b1, d, l);
    int len = 1;
    while (len < n + l) len <<= 1;
    ntt(a1, len, 1);
    ntt(d, len, 1);
    for (int i = 0; i < len; i++) d[i] = 1ll * a1[i] * d[i] % mod;
    ntt(d, len, -1);
    for (int i = l; i < len; i++) d[i] = 0;
    for (int i = 0; i < l >> 1; i++) swap(d[i], d[l - i - 1]);
    for (int i = 0; i < m; i++) r[i] = b[i];
    for (int i = 0; i < l; i++) b1[i] = d[i];
    for (int i = l; i < m; i++) b1[i] = 0;
    ntt(r, len, 1);
    ntt(b1, len, 1);
    for (int i = 0; i < len; i++) r[i] = 1ll * r[i] * b1[i] % mod;
    ntt(r, len, -1);
    for (int i = 0; i < n; i++) r[i] = (a[i] - r[i] + mod) % mod;
}
           

##求ln

##exp

這兩個不想寫了……留坑待填