Polynomial Related Operations & Calculations.

Table of Contents

  1. Polynomial Transform
    1. Fast Fourier Transformation
    2. Number-Theorem Transformation
    3. Arbitrary MOD Number-Theorem Transformation
  • Polynomial Operation
    1. Inversion
    2. Mod
    3. Log
    4. Exponent

    • Polynomial Transform

    Fast Fourier Transformation

    折半引理，消去引理，求和引理

    复数单位根

    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
    #include<iostream>
    #include<cstdio>
    #include<algorithm>
    #include<cmath>

    using std::acos;
    using std::cos;
    using std::sin;
    using std::swap;

    const double PI = acos(-1.0);
    const int maxn = 2500005;

    struct comp{
        double r,i;
        comp(double rr = 0.0, double ii = 0.0)
        {
            r = rr;
            i = ii;
        }
    }a[maxn], b[maxn], omega[maxn], omega_[maxn];
    comp operator + (comp a, comp b)
    {
        return comp(a.r + b.r, a.i + b.i);
    }
    comp operator - (comp a, comp b)
    {
        return comp(a.r - b.r, a.i - b.i);
    }
    comp operator * (comp a, comp b)
    {
        return comp(a.r * b.r - a.i * b.i, a.r * b.i + a.i * b.r);
    }

    void pre_omega(int x)
    {
        double round = 2.0*PI/x;
        for(int i = 0; i < x; ++i)
            omega[i] = comp(cos(i * round),sin(i * round)),
            omega_[i] = comp(cos(i * round),-1 * sin(i * round));//conj
    }
    void FFT(int n, comp *a, comp *w)
    {
        for(int i = 0, j = 0; i < n; ++i)
        {
            if(i > j) swap(a[i] , a[j]);
            for(int l = (n >> 1); (j ^= l) < l; l >>= 1); //remember it
        }
        for(int i = 2; i <= n; i <<= 1)
        {
            int m = (i >> 1);
            for(int j = 0; j < n; j += i)
                for(int k = 0; k < m; ++k)
                {
                    comp t = a[j + m + k] * w[n / i * k];
                    a[j + m + k] = a[j + k] - t;
                    a[j + k] = a[j + k] + t;
                }
                    
        }
    }

    int main()
    {
        int n, m, l = 1;
        scanf("%d%d",&n, &m);
        for(int i = 0; i <= n; ++i)
            scanf("%lf",&a[i].r);
        for(int i = 0; i <= m; ++i)
            scanf("%lf",&b[i].r);
        while(l <= n + m) l <<= 1;
        pre_omega(l);
        FFT(l, a, omega);
        FFT(l, b, omega);
        for(int i = 0; i <= l; ++i)
            a[i] = a[i] * b[i];
        FFT(l, a, omega_);
        for(int i = 0; i <= n + m; ++i) printf("%d ",(int)(a[i].r/l + 0.5));
        putchar('\n');
        return 0;
    }

    Number-Theorem Transformation

    在模意义下，质数原根具有和复数单位根同样的性质。

    与FFT流程大致相同，只是更换了原根，最后除以$ 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
    #include <iostream>
    #include <algorithm>
    #include <cstring>
    #include <cstdio>
    #include <cmath>
    #include <cctype>

    using namespace std;
    typedef long long ll;

    const ll mod = 998244353;
    const int MAX = 1e+5 + 5e+4 + 5;

    int n, m;
    char str[MAX];
    ll A[1 << 18], B[1 << 18], h[1 << 18];
    ll e[MAX], ine[MAX], inv3;

    inline ll quick_power(ll base, int index)
    {
    	ll ret = 1;
    	while(index)
    	{
    		if(index & 1) ret = (ret * base) % mod;
    		index >>= 1;
    		base = (base * base) % mod;
    	}
    	return ret;
    }
    inline void pre()
    {
    	ll inv3 = quick_power(3, mod - 2);
    	for(int i = 1; i <= n; i <<= 1) e[i] = quick_power(3, (mod - 1) / i);
    	for(int i = 1; i <= n; i <<= 1) ine[i] = quick_power(inv3, (mod - 1) / i);
    }
    inline void DFT(ll *a, int n, int t)
    {
    	for(int i = 0, j = 0; i < n; ++i)
    	{
    		if(i > j) swap(a[i], a[j]);
    		for(int l = (n >> 1); (j ^= l) < l; l >>= 1);
    	}
    	for(int i = 2; i <= n; i <<= 1)
    	{
    		int m = (i >> 1);
    		for(int j = 0; j < n; j += i)
    		{
    			ll w = 1, wn = e[i];
    			if(t == -1) wn = ine[i];
    			for(int k = 0; k < m; ++k)
    			{
    				ll z = a[j + m + k] * w % mod;
    				a[j + m + k] = (a[j + k] - z + mod) % mod;
    				a[j + k] = (a[j + k] + z) % mod;
    				w = (w * wn) % mod;
    			}
    		}
    	}
    }

    int main()
    {
    	scanf("%d", &n);
    	scanf("%s", str);
    	for(int i = 0; i < n; ++i) A[i] = str[n - i - 1] - '0';
    	scanf("%s", str);
    	for(int i = 0; i < n; ++i) B[i] = str[n - i - 1] - '0';
    	m = (n * 2), n = 1;
    	while(n <= m) n <<= 1;
    	pre();
    	DFT(A, n, 1); 
    	DFT(B, n, 1);
    	for(int i = 0; i < n; ++i) A[i] = A[i] * B[i] % mod;
    	DFT(A, n, -1);
    	ll inv = quick_power((ll)n, mod - 2);
    	for(int i = 0; i < n; ++i) A[i] = A[i] * inv % mod;
    	for(int i = 0; i < n; ++i) A[i + 1] += A[i]/10, A[i] %= 10;
    	int top = n - 1; while(top && !A[top]) top--;
    	for(int i = top; i >= 0; --i) putchar(A[i] + '0');
    	puts(""); 
    	return 0;
    }

    Arbitrary MOD Number-Theorem Transformation

    在任意模数下，很可能没有原根，但是如果要求精度很高，FFT又无法胜任。

    一个经典的做法是选取三个模数再用中国剩余定理合并。

    更高效的做法是利用上多项式的实部和虚部性质，把多项式系数拆分成前15位和后15位。做的运算更少。

    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 <cmath>
    #include <cstdio>

    using std::acos;
    using std::cos;
    using std::sin;
    using std::reverse;
    using std::swap;

    typedef long long ll;

    const int maxn = 262150;
    const int mask15 = 32767;
    const double pi = acos(-1);

    struct comp
    {
    	double r, i;
    	comp(double r_ = 0.0, double i_ = 0.0):r(r_),i(i_){}
    	comp operator + (const comp &rhs)
    	{
    		return comp(r + rhs.r, i + rhs.i);
    	}
    	comp operator - (const comp &rhs)
    	{
    		return comp(r - rhs.r, i - rhs.i);
    	}
    	comp operator * (const comp &rhs)
    	{
    		return comp(r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r);
    	}
    	
    };
    comp conj(const comp &c)
    {
    	return comp(c.r, -c.i);
    }

    int nn, mm, N, L;
    int mod;
    int rev[maxn];
    comp omega[maxn];

    inline void init()
    {
    	for(int i = 0; i < N; ++i)
    		omega[i] = comp(cos(2 * pi * i / N), sin(2 * pi * i / N)),
    		rev[i] = rev[i >> 1] >> 1 | ((i & 1) << (L - 1));
    }
    inline void FFT(int n, comp *a)
    {
    	for(int i = 0; i < n; ++i) if(i < rev[i]) swap(a[i], a[rev[i]]);
    	for(int i = 2, j = (n >> 1); i <= n; i <<= 1, j >>= 1)
    		for(int k = 0; k < n; k += i)
    		{
    			comp *l = a + k, *r = a + k + (i >> 1), *w = omega;
    			for(int p = 0; p < (i >> 1); ++p)
    			{
    				comp t = (*r) * (*w);
    				*r = (*l) - t, *l = (*l) + t;
    				++l; ++r; w += j;
    			} 
    		}
    }
    inline void solve(int *x, int *y, int *z)
    {
    	for(int i = 0; i < N; ++i) 
    		x[i] = (int)((ll)(x[i] + mod) % mod), y[i] = (int)((ll)(y[i] + mod) % mod);
    	static comp a[maxn], b[maxn];
    	static comp dft1[maxn], dft2[maxn], dft3[maxn], dft4[maxn];

    	for(int i = 0; i < N; ++i) 
    		a[i] = comp(x[i] & mask15, x[i] >> 15),
    		b[i] = comp(y[i] & mask15, y[i] >> 15);
    	FFT(N, a); FFT(N, b);
    	for(int i = 0; i < N; ++i)
    	{
    		int j = (N - i) & (N - 1);
    		static comp da, db, dc, dd;
    		da = (a[i] + conj(a[j])) * comp(0.5, 0.0);
    		db = (a[i] - conj(a[j])) * comp(0.0, -0.5);
    		dc = (b[i] + conj(b[j])) * comp(0.5, 0.0);
    		dd = (b[i] - conj(b[j])) * comp(0.0, -0.5);
    		dft1[j] = da * dc;
    		dft2[j] = da * dd;
    		dft3[j] = db * dc;
    		dft4[j] = db * dd;
    	}
    	for(int i = 0; i < N; ++i) a[i] = dft1[i] + dft2[i] * comp(0, 1);
    	for(int i = 0; i < N; ++i) b[i] = dft3[i] + dft4[i] * comp(0, 1);
    	FFT(N, a); FFT(N, b);
    	for(int i = 0; i < N; ++i)
    	{
    		int da = (ll)(a[i].r / N + 0.5) % mod;
    		int db = (ll)(a[i].i / N + 0.5) % mod;
    		int dc = (ll)(b[i].r / N + 0.5) % mod;
    		int dd = (ll)(b[i].i / N + 0.5) % mod;
    		z[i] = (da + ((ll)(db + dc) << 15) + ((ll)dd << 30)) % mod;
    	}
    }

    int main()
    {
    	scanf("%d%d%d", &nn, &mm, &mod); nn++; mm++;
    	static int A[maxn], B[maxn], C[maxn];
    	for(int i = 0; i < nn; ++i) 
    		scanf("%d", &A[i]);
    	for(int i = 0; i < mm; ++i)
    		scanf("%d", &B[i]);
    	L = 0;
    	for(; (1 << L) < nn + mm - 1; ++L);
    	N = (1 << L);
    	init();
    	solve(A, B, C);
    	for(int i = 0; i < nn + mm - 1; ++i) 
    		C[i] = (C[i] + mod) % mod, printf("%d ", C[i]);
    	return 0;
    }

    Polynomial Operation

    Inversion

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    void Inversion(int *arr, int *res, int len) {
      if (len == 1) {
        res[0] = quick_power(arr[0], mod - 2);
        return;
      }
      Inversion(arr, res, len >> 1);
      int now = len << 1;
      for (int i = 0; i < len; ++i) tmp[i] = arr[i];
      for (int i = len; i < now; ++i) tmp[i] = 0;
      DFT(tmp, 1, now); DFT(res, 1, now);
      for (int i = 0; i < now; ++i)
        tmp[i] = 1ll * res[i] * (2ll - 1ll * tmp[i] * res[i] % mod + mod) % mod;
      DFT(tmp, -1, now);
      for (int i = 0; i < len; ++i) res[i] = tmp[i];
      for (int i = len; i < now; ++i) res[i] = 0;
    }

    Mod

    先咕着.

    Log

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    void Logarithmic(int *arr, int *res, int len) {
      Derivative(arr, tmp, len);
      Inversion(arr, tmp2, len);
      memset(tmp, 0, sizeof tmp);
      memset(tmp2, 0, sizeof tmp);
      DFT(tmp, 1, len << 1); DFT(tmp2, 1, len << 1);
      for (int i = 0; i < (len << 1); ++i)
        tmp[i] = 1ll * tmp[i] * tmp2[i] % mod;
      DFT(tmp, -1, len << 1);
      Intergral(tmp, res, len);
    }

    Exponent

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    void Exponent(int *arr, int *res, int len) {
      if (len == 1) {
        res[0] = 1;
        return;
      }
      Exponent(arr, res, len >> 1);
      for (int i = 0; i < len; ++i) tmp3[i] = res[i];
      Logarithmic(tmp3, tmp4, len);
      for (int i = 0; i < len; ++i)
        tmp4[i] = (arr[i] - tmp4[i] + mod) % mod;
      tmp4[0] = (tmp4[0] + 1) % mod;
      DFT(tmp3, 1, len << 1); DFT(tmp4, 1, len << 1);
      for (int i = 0; i < (len << 1); ++i) tmp3[i] = 1ll * tmp3[i] * tmp4[i] % mod;
      DFT(tmp3, -1, (len << 1));
      for (int i = 0; i < len; ++i) res[i] = tmp3[i];
      memset(tmp3, 0, sizeof tmp);
      memset(tmp4, 0, sizeof tmp2);
    }