Loading [MathJax]/jax/output/CommonHTML/jax.js
My Algorithm : kopricky アルゴリズムライブラリ

kopricky アルゴリズムライブラリ

Polynomial Logarithm

コードについての説明

体 K 上の多項式 f (deg(f)=n, f(0)=1) に対応する形式的べき級数 F が与えられたときに 形式的べき級数 loge(F) の xi(0≤i<r) の係数を返す.
以下の実装は K=Z/pZ (p は素数) の場合の実装で特に素数 p について p=k∗2l+1(k>0,n+1≤2l) を満たすこと(数論変換に乗ること)を仮定している.
正直形式的べき級数をあまり良く理解していないが, ちゃんと理解することへのモチベーションが湧いていないので以下雰囲気で説明を書く.
FPS についての記事 によると loge(F) は loge(F)=∫F′F の関係を用いることで簡単に計算できる(ここでの = は形式的べき級数の意味での等号).
無限級数の積分を見ると, 一様収束するなら項別積分可能だけど...とかいう気持ちになるが形式的べき級数はそういう概念ではないっぽい(？)ので 可積分性とかは特に考えないことにする.

(関数)
polynomial_log(a,r): 多項式 f の係数 a(a[0]=1) を与えたとき, それに対応する形式的べき級数 F とする. このとき loge(F) の xi(0≤i<r) の係数を返す.

時間計算量: O(n+rlog2r)

コード

  1. #define MOD 998244353
  2. #define root 3
  3.  
  4. unsigned int add(const unsigned int x, const unsigned int y)
  5. {
  6.     return (x + y < MOD) ? x + y : x + y - MOD;
  7. }
  8.  
  9. unsigned int sub(const unsigned int x, const unsigned int y)
  10. {
  11.     return (x >= y) ? (x - y) : (MOD - y + x);
  12. }
  13.  
  14. unsigned int mul(const unsigned int x, const unsigned int y)
  15. {
  16.     return (unsigned long long)x * y % MOD;
  17. }
  18.  
  19. unsigned int mod_pow(unsigned int x, unsigned int n)
  20. {
  21.     unsigned int res = 1;
  22.     while(n > 0){
  23.         if(n & 1){ res = mul(res, x); }
  24.         x = mul(x, x);
  25.         n >>= 1;
  26.     }
  27.     return res;
  28. }
  29.  
  30. unsigned int inverse(const unsigned int x)
  31. {
  32.     return mod_pow(x, MOD - 2);
  33. }
  34.  
  35. void ntt(vector<int>& a, const bool rev = false)
  36. {
  37.     unsigned int i, j, k, l, p, q, r, s;
  38.     const unsigned int size = a.size();
  39.     if(size == 1) return;
  40.     vector<int> b(size);
  41.     r = rev ? (MOD - 1 - (MOD - 1) / size) : (MOD - 1) / size;
  42.     s = mod_pow(root, r);
  43.     vector<unsigned int> kp(size / 2 + 1, 1);
  44.     for(i = 0; i < size / 2; ++i) kp[i + 1] = mul(kp[i], s);
  45.     for(i = 1, l = size / 2; i < size; i <<= 1, l >>= 1){
  46.         for(j = 0, r = 0; j < l; ++j, r += i){
  47.             for(k = 0, s = kp[i * j]; k < i; ++k){
  48.                 p = a[k + r], q = a[k + r + size / 2];
  49.                 b[k + 2 * r] = add(p, q);
  50.                 b[k + 2 * r + i] = mul(sub(p, q), s);
  51.             }
  52.         }
  53.         swap(a, b);
  54.     }
  55.     if(rev){
  56.         s = inverse(size);
  57.         for(i = 0; i < size; i++){ a[i] = mul(a[i], s); }
  58.     }
  59. }
  60.  
  61. vector<int> convolute(const vector<int>& a, const vector<int>& b, int asize, int bsize, int _size)
  62. {
  63.     const int size = asize + bsize - 1;
  64.     int t = 1;
  65.     while (t < size){ t <<= 1; }
  66.     vector<int> A(t, 0), B(t, 0);
  67.     for(int i = 0; i < asize; i++){ A[i] = (a[i] < MOD) ? a[i] : (a[i] % MOD); }
  68.     for(int i = 0; i < bsize; i++){ B[i] = (b[i] < MOD) ? b[i] : (b[i] % MOD); }
  69.     ntt(A), ntt(B);
  70.     for(int i = 0; i < t; i++){ A[i] = mul(A[i], B[i]); }
  71.     ntt(A, true);
  72.     A.resize(_size);
  73.     return A;
  74. }
  75.  
  76. vector<int> polynomial_inverse(const vector<int>& a, int r)
  77. {
  78.     assert(a[0] != 0);
  79.     vector<int> h = {(int)inverse(a[0])};
  80.     int t = 1;
  81.     for(int i = 0; t < r; ++i){
  82.         t <<= 1;
  83.         vector<int> res = convolute(a, convolute(h, h, t / 2, t / 2, t), min((int)a.size(), t), t, t);
  84.         for(int j = 0; j < t; ++j){
  85.             res[j] = sub(0, res[j]);
  86.             if(j < t / 2) res[j] = add(res[j], mul(2, h[j]));
  87.         }
  88.         swap(h, res);
  89.     }
  90.     h.resize(r);
  91.     return h;
  92. }
  93.  
  94. vector<int> differentiate(const vector<int>& a)
  95. {
  96.     const int n = (int)a.size();
  97.     vector<int> res(n - 1);
  98.     for(int i = 0; i < n - 1; ++i){
  99.         res[i] = mul(a[i + 1], i + 1);
  100.     }
  101.     return res;
  102. }
  103.  
  104. vector<int> integral(const vector<int>& a)
  105. {
  106.     const int n = (int)a.size();
  107.     vector<int> res(n + 1, 0);
  108.     for(int i = 0; i < n; ++i){
  109.         res[i + 1] = mul(a[i], inverse(i + 1));
  110.     }
  111.     return res;
  112. }
  113.  
  114. vector<int> polynomial_log(const vector<int>& a, int r)
  115. {
  116.     assert(a[0] == 1);
  117.     if(r == 0) return {0};
  118.     const int n = (int)a.size();
  119.     const vector<int> num = differentiate(a);
  120.     const vector<int> den = polynomial_inverse(a, r - 1);
  121.     const vector<int> res = convolute(num, den, min(n, r) - 1, r - 1, r - 1);
  122.     return integral(res);
  123. }

verify 用の問題

yosupo さんの library checker : Log of Formal Power Series 提出コード