多項式相關算法集成

時間 2019-11-11

標籤多項式相關算法集成简体版

原文原文鏈接

原理

快速數論變換(FNT)

原理同FFT同樣,用與弧度具備一樣性質的一組值替代弧度進行計算,具體性質以下:html

FFT 中能在 O(nlog⁡n) 時間內變換用到了單位根 \(\omega\)的什麼性質:ios

\(\omega_n^n = 1\)c++
\(\omega_{n}^{0}\), \(\omega_{n}^{1}\), ⋯, \(\omega_{n}^{n-1}\)是互不相同的，這樣帶入計算出來的點值才能夠用來還原出係數算法
\(\omega_{n}^{2}=\omega_{n/2}^{1}\), \(\omega_{n}^{n/2+k}=-\omega_{n}^{k}\)，這使得在按照指數奇偶分類以後可以把帶入的值也減半使得問題規模可以減半數組
\[ \sum_{k=0}^{n-1}(\omega_{n}^{j-i})^{k}=\begin{cases} 0,& \text{若i != j}\\ 1,& \text{若i = j} \end{cases} \]spa

這點保證了可以使用相同的方法進行逆變換獲得係數表示code

參考:orm

多項式逆元

參考:

多項式求逆元

多項式除法及求模

參考:

多項式除法及求模

多項式多點求值與拉格朗日插值

參考:

實現

展開查看

/*
 * 多項式相關算法模板集合
 */
#include 
   
   
   

  
   
  using namespace std; typedef long long ll; // 2281701377是一個原根爲3的指數,平方恰好不會爆long long // 另一個經常使用的998244353 = 119ll * (1 << 23) + 1 // 1004535809=479⋅221+1 加起來恰好不會爆 int 也不錯 const ll mod_v = 17ll * (1 << 27) + 1; const int MAXL = 18, N = 1 << MAXL, M = 1 << MAXL; // MAXL最大取mov_v的2次冪 struct polynomial {  ll coef[N]; // size設爲多項式係數個數的兩倍  ll a[N], b[N], d[N], r[N];  ll xcoor[N], ycoor[N], v[N], mcoef[N]; // size設爲點個數的兩倍  vector 
 
    
    
      > poly_divisor;   static ll mod_pow(ll x,ll n,ll m) {  ll ans = 1;  while(n > 0){  if(n & 1)  ans = ans*x%m;  x = x*x%m;  n >>= 1;  }  return ans;  }  struct FastNumberTheoreticTransform { ll omega[N], omegaInverse[N]; int range; // 初始化頻率 void init (const int& n) { range = n; ll base = mod_pow(3, (mod_v - 1) / n, mod_v); ll inv_base = mod_pow(base, mod_v - 2, mod_v); omega[0] = omegaInverse[0] = 1; for (int i = 1; i < n; ++i) { omega[i] = omega[i - 1] * base % mod_v; omegaInverse[i] = omegaInverse[i - 1] * inv_base % mod_v; } } // Cooley-Tukey算法:O(n*logn) void transform (ll *a, const ll *omega, const int &n) { for (int i = 0, j = 0; i < n; ++i) { if (i > j) std::swap (a[i], a[j]); for(int l = n >> 1; ( j ^= l ) < l; l >>= 1); } for (int l = 2; l <= n; l <<= 1) { int m = l / 2; for (ll *p = a; p != a + n; p += l) { for (int i = 0; i < m; ++i) { ll t = omega[range / l * i] * p[m + i] % mod_v; p[m + i] = (p[i] - t + mod_v) % mod_v; p[i] = (p[i] + t) % mod_v; } } } } // 時域轉頻域 void dft (ll *a, const int& n) { transform(a, omega, n); } // 頻域轉時域 void idft (ll *a, const int& n) { transform(a, omegaInverse, n); for (int i = 0; i < n; ++i) a[i] = a[i] * mod_pow(n, mod_v - 2, mod_v) % mod_v; } } fnt ; // 與模比較轉換值 ll mod_trans(ll v) { return abs(v) <= mod_v / 2 ? v : (v < 0 ? v + mod_v : v - mod_v); } // 二分求\prod{x-xi}的全部二分子多項式的係數 void binary_subpoly(int l, int r, int idx) { if (l == r - 1) { poly_divisor[idx].push_back(-xcoor[l]); poly_divisor[idx].push_back(1); } else { int lidx = (idx << 1) + 1, ridx = lidx + 1; binary_subpoly(l, (l + r) / 2, lidx); binary_subpoly((l + r) / 2, r, ridx); int t = poly_divisor[lidx].size() + poly_divisor[ridx].size() - 1; int p = 1; while(p < t) p <<= 1; copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a); fill(a + poly_divisor[lidx].size(), a + p, 0); copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b); fill(b + poly_divisor[ridx].size(), b + p, 0); fnt.dft(a, p); fnt.dft(b, p); for (int i = 0; i < p; i++) a[i] *= b[i]; fnt.idft(a, p); for (int i = 0; i < t; i++) poly_divisor[idx].push_back(mod_trans(a[i])); } } // 模x^deg,a爲要求逆元的多項式係數,結果存放在b[0~deg]中 // T(deg) = T(deg/2) + deg*log(deg)，複雜度O(deg*log(deg)) void polynomial_inverse(int deg, ll* a, ll* b, ll* tmp) { if(deg == 1) { b[0] = mod_pow(a[0], mod_v - 2, mod_v); } else { polynomial_inverse((deg + 1) >> 1, a, b, tmp); int p = 1; while(p < (deg << 1) - 1) p <<= 1; copy(a, a + deg, tmp); fill(tmp + deg, tmp + p, 0); fill(b + ((deg + 1) >> 1), b + p, 0); //fnt.init(p); fnt.dft(tmp, p); fnt.dft(b, p); for(int i = 0; i != p; ++i) { b[i] = (2 - tmp[i] * b[i] % mod_v) * b[i] % mod_v; if(b[i] < 0) b[i] += mod_v; } fnt.idft(b, p); fill(b + deg, b + p, 0); } } // A = D*B + R,A爲n項n-1次冪,B爲m項m-1次冪,D爲n-m+1項n-m次冪,R爲m-1項m-2次冪 // 要求a,b中係數以低次到屢次順序排列 // n >= m，複雜度O(n*logn)；n < m，複雜度O(n) int polynomial_division(int n, int m, ll *A, ll *B, ll *D, ll *R) { if (n < m) { copy(A, A + n, R); return n; } else { static ll A0[N], B0[N], tmp[N]; //數組太大會爆棧,添加到全局區 int p = 1, t = n - m + 1; while(p < (t << 1) - 1) p <<= 1; fill(A0, A0 + p, 0); reverse_copy(B, B + m, A0); polynomial_inverse(t, A0, B0, tmp); fill(B0 + t, B0 + p, 0); fnt.dft(B0, p); reverse_copy(A, A + n, A0); fill(A0 + t, A0 + p, 0); fnt.dft(A0, p); for(int i = 0; i != p; ++i) A0[i] = A0[i] * B0[i] % mod_v; fnt.idft(A0, p); reverse(A0, A0 + t); copy(A0, A0 + t, D); for(p = 1; p < n; p <<= 1); fill(A0 + t, A0 + p, 0); fnt.dft(A0, p); copy(B, B + m, B0); fill(B0 + m, B0 + p, 0); fnt.dft(B0, p); for(int i = 0; i != p; ++i) A0[i] = A0[i] * B0[i] % mod_v; fnt.idft(A0, p); for(int i = 0; i != m - 1; ++i) R[i] = ((A[i] - A0[i]) % mod_v + mod_v) % mod_v; //fill(R + m - 1, R + p, 0); return m - 1; } } // 多項式的點值計算 // l和r爲存儲要求的點數組的左右邊界(左閉右開), idx爲除數多項式索引(初始０) // polycoef爲用於計算點的多項式係數，num爲其係數個數 // 設多項式項數x=num，點數y=r-l，n=max(x,y)，複雜度O(n(logn)^2) void polynomial_calculator(int l, int r, int idx, int num, ll *polycoef) { int mid = (l + r) / 2; int lidx = (idx << 1) + 1, ridx = lidx + 1; int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size(); ll *lmod_poly = new ll[lsize - 1], *rmod_poly = new ll[rsize - 1]; copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a); copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b); int lplen = polynomial_division(num, lsize, polycoef, a, d, lmod_poly); int rplen = polynomial_division(num, rsize, polycoef, b, d, rmod_poly); if (l == mid - 1) { v[l] = mod_trans(lmod_poly[0]); } else { polynomial_calculator(l, (l + r) / 2, lidx, lplen, lmod_poly); } if (r == mid + 1) { v[(l + r) / 2] = mod_trans(rmod_poly[0]); } else { polynomial_calculator((l + r) / 2, r, ridx, rplen, rmod_poly); } delete []lmod_poly; delete []rmod_poly; } // 拉格朗日插值：二分＋快速數論變換 // l和r爲存儲要求的點數組的左右邊界(左閉右開), idx爲由點二分構造出多項式的索引(初始０) // polycoef爲插值獲得的多項式結果 // 設點個數爲n，複雜度O(n*(logn)^2)，結果polycoef爲n-1次多項式 void polynomial_interpolate(int l, int r, int idx, ll *polycoef) { if (l == r - 1) { polycoef[0] = ycoor[l] * mod_pow(v[l], mod_v - 2, mod_v) % mod_v; } else { int mid = (l + r) >> 1; int lidx = (idx << 1) + 1, ridx = lidx + 1; int sz = poly_divisor[idx].size() - 1; int lsize = poly_divisor[lidx].size(), rsize = poly_divisor[ridx].size(); int p = 1; while (p < sz) p <<= 1; ll *leftpoly = new ll[p], *rightpoly = new ll[p]; polynomial_interpolate(l, mid, lidx, leftpoly); polynomial_interpolate(mid, r, ridx, rightpoly); copy(poly_divisor[lidx].begin(), poly_divisor[lidx].end(), a); copy(poly_divisor[ridx].begin(), poly_divisor[ridx].end(), b); fill(leftpoly + lsize - 1, leftpoly + p, 0); fill(rightpoly + rsize - 1, rightpoly + p, 0); fill(a + lsize, a + p, 0); fill(b + rsize, b + p, 0); fnt.dft(leftpoly, p); fnt.dft(b, p); for (int i = 0; i < p; i++) leftpoly[i] = leftpoly[i] * b[i] % mod_v; fnt.idft(leftpoly, p); fnt.dft(rightpoly, p); fnt.dft(a, p); for (int i = 0; i < p; i++) rightpoly[i] = rightpoly[i] * a[i] % mod_v; fnt.idft(rightpoly, p); for (int i = 0; i < sz; i++) polycoef[i] = mod_trans((leftpoly[i] + rightpoly[i]) % mod_v); delete []leftpoly; delete []rightpoly; } } // 初始化點個數到二分子多項式個數 // 調用polynomial_calculator和polynomial_interpolate前調用 void init(int vnum) { int vnum2 = 1; while (vnum2 < vnum) vnum2 <<= 1; for (int i = 0; i < 2 * vnum2 - 1; i++) poly_divisor.push_back(vector 
     
       ()); } } poly;

應用

多項式快速乘

展開代碼

int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, m;
    cin >> n >> m;
    for (int i = 0; i < n; i++) cin >> poly.a[i];
    for (int i = 0; i < m; i++) cin >> poly.b[i];
    cout << "a*b以後的真實係數:" << endl;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < m; j++)
            poly.r[i + j] += poly.a[i] * poly.b[j];
    for (int i = 0; i < n + m - 1; i++)
        cout << poly.r[i] << " ";
    cout << endl;
    cout << "利用fft計算得出的係數:" << endl;
    int p = 1;
    while (p < n + m - 1) p <<= 1;  // 只要p大於多項式結果中的最大次冪便可
    poly.fnt.dft(poly.a, p);
    poly.fnt.dft(poly.b, p);
    for (int i = 0; i < p; i++) {
        poly.d[i] = poly.a[i] * poly.b[i] % mod_v;
    }
    poly.fnt.idft(poly.d, p);
    for (int i = 0; i < n + m - 1; i++) cout << poly.d[i] << " ";
    return 0;
}

多項式逆元

展開代碼

int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n;
    cin >> n;
    for(int i = 0; i != n; ++i)
        cin >> poly.a[i];
    int m;
    cin >> m;  // 輸入模x^m
    poly.polynomial_inverse(m, poly.a, poly.b, poly.d);
    cout << "inverse: " << endl;
    for(int i = 0; i != m; ++i)
        cout << (poly.b[i] + mod_v) % mod_v << " ";
    cout << endl;
    cout << "a*b相乘取模驗證取模後的前m項係數:" << endl;
    memset(poly.d, 0, sizeof(poly.d));
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            poly.d[i + j] = (poly.d[i + j] + poly.a[i] * poly.b[j] % mod_v)% mod_v;
        }   
    }   
    cout << m << endl;
    for (int i = 0; i < m; i++) {
        cout << poly.d[i] << " ";
    }   
    return 0;
}

輸入數據獲得輸出後的結果以下圖

多項式除法以及取模

展開代碼

int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, m;
    cin >> n >> m;
    for(int i = 0; i != n; ++i)
    cin >> poly.a[i]; // 0次冪係數開始輸入,缺失的冪係數輸入0
    for (int i = 0; i < m; i++)
    cin >> poly.b[i]; // 0次冪係數開始輸入
    int rlen = poly.polynomial_division(n, m, poly.a, poly.b, poly.d, poly.r);
    cout << "輸出商多項式:" << endl;
    for (int i = 0; i < n - m + 1; i++) 
    cout << poly.d[i] << " ";
    cout << endl;
    cout << "輸出餘數多項式:" << endl;
    for (int i = 0; i < rlen; i++)
    cout << poly.r[i] << " ";
    return 0;
}

多項式多點求值

展開代碼

int main() {
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);
    int n, vnum; 
    // 輸入要計算的點
    cin >> n;
    for (int i = 0; i < n; i++) {
        cin >> poly.coef[i];
    }
    cin >> vnum;
    for (int i = 0; i < vnum; i++) {
        cin >> poly.xcoor[i];        
    }
    // 二分求子多項式                                                                               
    poly.init(vnum);
    poly.binary_subpoly(0, vnum / 2, 1);
    poly.binary_subpoly(vnum / 2, vnum, 2); 
    // 將輸入點代入插值獲得的多項式中進行驗證，輸出計算獲得的結果
    cout << "計算結果:" << endl;
    poly.polynomial_calculator(0, vnum, 0, n, poly.coef);
    for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
    cout << endl;
}

輸入數據後獲得結果以下

拉格朗日快速插值計算

展開代碼

int main() {
    // 多項式係數默認低次到高次排列
    ios_base::sync_with_stdio(false);
    poly.fnt.init(1 << MAXL);   
    // 多項式插值
    int vnum; 
    // 輸入要計算的點
    cin >> vnum;
    for (int i = 0; i < vnum; i++) {
        cin >> poly.xcoor[i] >> poly.ycoor[i];        
    }
    // 二分求子多項式
    poly.init(vnum);
    poly.binary_subpoly(0, vnum, 0); 
    for (unsigned int i = 1; i < poly.poly_divisor[0].size(); i++) {
        poly.mcoef[i - 1] = poly.poly_divisor[0][i] * i % mod_v;
    }
    // 遍歷i計算全部\sum_{j!=i}{xi-xj} 
    poly.polynomial_calculator(0, vnum, 0, poly.poly_divisor[0].size() - 1, poly.mcoef);
    // 拉格朗日插值計算多項式
    poly.polynomial_interpolate(0, vnum, 0, poly.coef);
    // 輸出插值獲得的多項式係數 
    for (int i = 0; i < vnum; i++) {
        cout << poly.coef[i] << " ";
    }
    cout << endl;
    // 將輸入點代入插值獲得的多項式中進行驗證，輸出計算獲得的結果
    poly.polynomial_calculator(0, vnum, 0, vnum, poly.coef);
    for (int i = 0; i < vnum; i++) cout << poly.v[i] << " ";
    cout << endl;
    return 0;
}

輸入數據獲得的結果以下