fork download 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3.  
  4. // ===================================================================
  5. // This file contains a collection of Fast Fourier Transform (FFT) and
  6. // Number Theoretic Transform (NTT) algorithms. Each function is ready
  7. // to be used as a "black box".
  8. // Read the comments above each one to understand:
  9. //   - What it solves
  10. //   - What input it expects
  11. //   - What it returns
  12. //   - Time complexity
  13. //   - Important constraints / assumptions
  14. // ===================================================================
  15.  
  16. // ===================================================================
  17. // 1) Core NTT (Number Theoretic Transform) Implementation
  18. //    NTT is the modular arithmetic version of FFT. It computes the
  19. //    exact convolution of two integer sequences modulo a prime number.
  20. //    This is the foundation for almost all functions in this file.
  21. //    Think of it as: given two polynomials, NTT multiplies them very fast.
  22. // ===================================================================
  23.  
  24. // NTT-friendly primes and their primitive roots.
  25. // A "primitive root" is a special number that generates all non-zero
  26. // elements of the finite field when raised to different powers.
  27. // The modulus must be of the form: mod = c * 2^k + 1, where k is large
  28. // enough to support the transform length.
  29. const int MOD = 998244353;        // = 119 * 2^23 + 1, primitive root = 3
  30. const int MOD2 = 1004535809;      // = 479 * 2^21 + 1, primitive root = 3
  31. const int MOD3 = 469762049;       // = 7 * 2^26 + 1, primitive root = 3
  32. const int PRIMITIVE_ROOT = 3;
  33.  
  34. // 1.1) Modular exponentiation (fast power).
  35. //      This is a helper function used internally by NTT.
  36. //      Parameters:
  37. //        - a: base (long long)
  38. //        - e: exponent (long long)
  39. //        - mod: modulus (long long)
  40. //      Returns:
  41. //        - a^e % mod
  42. //      Time complexity: O(log e)
  43. //      Constraint: mod > 0.
  44. long long modpow(long long a, long long e, long long mod) {
  45.     long long r = 1;
  46.     while (e) {
  47.         if (e & 1) r = (r * a) % mod;
  48.         a = (a * a) % mod;
  49.         e >>= 1;
  50.     }
  51.     return r;
  52. }
  53.  
  54. // 1.2) Generic NTT (works with any NTT‑friendly modulus and primitive root).
  55. //      This is the core transform function. It converts a polynomial from
  56. //      coefficient form to point-value form (or vice versa) using the
  57. //      NTT algorithm.
  58. //      Think of it as a fast way to evaluate a polynomial at many points.
  59. //      Parameters:
  60. //        - a: vector of integers (the polynomial coefficients). This vector
  61. //             is modified in-place.
  62. //        - invert: boolean. If false, performs forward NTT.
  63. //                  If true, performs inverse NTT.
  64. //        - mod: the prime modulus (must be NTT‑friendly).
  65. //        - root: the primitive root modulo 'mod'.
  66. //      Returns:
  67. //        - Nothing (the result is stored in the input vector 'a').
  68. //      Time complexity: O(n log n), where n = a.size()
  69. //      Important constraints:
  70. //        - The length of 'a' (n) MUST be a power of two.
  71. //        - All coefficients must be in the range [0, mod-1].
  72. //      Note: This function uses the iterative Cooley-Tukey algorithm.
  73. void ntt_generic(vector<int>& a, bool invert, int mod, int root) {
  74.     int n = (int)a.size();
  75.  
  76.     // Bit-reversal permutation.
  77.     for (int i = 1, j = 0; i < n; i++) {
  78.         int bit = n >> 1;
  79.         for (; j & bit; bit >>= 1) j ^= bit;
  80.         j ^= bit;
  81.         if (i < j) swap(a[i], a[j]);
  82.     }
  83.  
  84.     for (int len = 2; len <= n; len <<= 1) {
  85.         int wlen = modpow(root, (mod - 1) / len, mod);
  86.         if (invert) wlen = modpow(wlen, mod - 2, mod);
  87.  
  88.         for (int i = 0; i < n; i += len) {
  89.             long long w = 1;
  90.             for (int j = 0; j < len / 2; j++) {
  91.                 int u = a[i + j];
  92.                 int v = (int)(a[i + j + len / 2] * w % mod);
  93.  
  94.                 a[i + j] = u + v;
  95.                 if (a[i + j] >= mod) a[i + j] -= mod;
  96.  
  97.                 a[i + j + len / 2] = u - v;
  98.                 if (a[i + j + len / 2] < 0) a[i + j + len / 2] += mod;
  99.  
  100.                 w = w * wlen % mod;
  101.             }
  102.         }
  103.     }
  104.  
  105.     if (invert) {
  106.         int n_inv = modpow(n, mod - 2, mod);
  107.         for (int &x : a) x = (int)((long long)x * n_inv % mod);
  108.     }
  109. }
  110.  
  111. // 1.3) NTT Convolution with a specified modulus and primitive root.
  112. //      Multiplies two polynomials modulo a given NTT-friendly prime.
  113. //      Parameters:
  114. //        - a, b: vectors of coefficients (values in [0, mod-1]).
  115. //        - mod: the NTT-friendly modulus.
  116. //        - root: the primitive root modulo 'mod'.
  117. //      Returns:
  118. //        - A vector of integers representing (a * b) modulo 'mod'.
  119. //      Time complexity: O(n log n) where n is the padded size.
  120. vector<int> convolution_mod(const vector<int>& a, const vector<int>& b, int mod, int root) {
  121.     int n = (int)a.size(), m = (int)b.size();
  122.     if (!n || !m) return {};
  123.  
  124.     int sz = 1;
  125.     while (sz < n + m - 1) sz <<= 1;
  126.  
  127.     vector<int> fa(a.begin(), a.end()), fb(b.begin(), b.end());
  128.     fa.resize(sz);
  129.     fb.resize(sz);
  130.  
  131.     ntt_generic(fa, false, mod, root);
  132.     ntt_generic(fb, false, mod, root);
  133.  
  134.     for (int i = 0; i < sz; i++) {
  135.         fa[i] = (int)((long long)fa[i] * fb[i] % mod);
  136.     }
  137.  
  138.     ntt_generic(fa, true, mod, root);
  139.     fa.resize(n + m - 1);
  140.     return fa;
  141. }
  142.  
  143. // 1.4) Default NTT Convolution (using the global MOD = 998244353).
  144. //      This is the simplest version; use it when your modulus is 998244353.
  145. vector<int> convolution(const vector<int>& a, const vector<int>& b) {
  146.     return convolution_mod(a, b, MOD, PRIMITIVE_ROOT);
  147. }
  148.  
  149. // 1.5) Arbitrary Modulus Convolution (using CRT with 3 NTT-friendly primes).
  150. //      This function multiplies two polynomials modulo an arbitrary integer 'mod'.
  151. //      It uses three NTT-friendly primes and combines the results via Garner's algorithm.
  152. //      Parameters:
  153. //        - a, b: vectors of coefficients (values in [0, mod-1]).
  154. //        - mod: the target modulus (can be any positive integer).
  155. //      Returns:
  156. //        - A vector of integers representing (a * b) modulo 'mod'.
  157. //      Time complexity: O(n log n) with a constant factor (~3x slower).
  158. //      Important constraints:
  159. //        - The product of the three primes is > 1e27, which is large enough for
  160. //          most practical ranges to avoid ambiguity.
  161. vector<int> convolutionArbitraryMod(const vector<int>& a, const vector<int>& b, int mod) {
  162.     const int m1 = 998244353, m2 = 1004535809, m3 = 469762049;
  163.     const int r1 = 3, r2 = 3, r3 = 3;
  164.  
  165.     auto c1 = convolution_mod(a, b, m1, r1);
  166.     auto c2 = convolution_mod(a, b, m2, r2);
  167.     auto c3 = convolution_mod(a, b, m3, r3);
  168.  
  169.     int n = (int)c1.size();
  170.     vector<int> res(n);
  171.  
  172.     // Precompute inverses for Garner's algorithm.
  173.     static const long long inv_m1_m2 = modpow(m1, m2 - 2, m2);
  174.     static const long long inv_m1m2_m3 = modpow((long long)m1 * m2 % m3, m3 - 2, m3);
  175.  
  176.     for (int i = 0; i < n; i++) {
  177.         long long x1 = c1[i];
  178.         long long t1 = ((c2[i] - x1) % m2 + m2) % m2;
  179.         t1 = t1 * inv_m1_m2 % m2;
  180.         long long x2 = x1 + (long long)m1 * t1;
  181.  
  182.         long long t2 = ((c3[i] - x2) % m3 + m3) % m3;
  183.         t2 = t2 * inv_m1m2_m3 % m3;
  184.         long long x3 = x2 + (long long)m1 * m2 % mod * t2 % mod;
  185.  
  186.         res[i] = (int)(x3 % mod);
  187.     }
  188.  
  189.     return res;
  190. }
  191.  
  192. // ===================================================================
  193. // 2) FFT (Fast Fourier Transform) with Complex Numbers
  194. //    FFT uses complex numbers and works for any real/integer input.
  195. //    It is more flexible than NTT (no modulus restrictions) but can
  196. //    have precision issues with very large numbers.
  197. // ===================================================================
  198.  
  199. using cd = complex<double>;
  200. const double PI = acos(-1);
  201.  
  202. // 2.1) Iterative FFT.
  203. //      This is the complex-number version of the transform.
  204. //      Parameters:
  205. //        - a: vector of complex numbers. Modified in-place.
  206. //        - invert: boolean. If false, forward transform. If true, inverse.
  207. //      Returns:
  208. //        - Nothing (result stored in 'a').
  209. //      Time complexity: O(n log n), where n = a.size()
  210. //      Important constraints:
  211. //        - The length of 'a' MUST be a power of two.
  212. //      Note: Uses complex numbers and may have precision errors.
  213. void fft(vector<cd>& a, bool invert) {
  214.     int n = a.size();
  215.  
  216.     // Bit-reversal permutation.
  217.     for (int i = 1, j = 0; i < n; i++) {
  218.         int bit = n >> 1;
  219.         for (; j & bit; bit >>= 1) j ^= bit;
  220.         j ^= bit;
  221.         if (i < j) swap(a[i], a[j]);
  222.     }
  223.  
  224.     for (int len = 2; len <= n; len <<= 1) {
  225.         double ang = 2 * PI / len * (invert ? -1 : 1);
  226.         cd wlen(cos(ang), sin(ang));
  227.  
  228.         for (int i = 0; i < n; i += len) {
  229.             cd w(1);
  230.             for (int j = 0; j < len / 2; j++) {
  231.                 cd u = a[i + j];
  232.                 cd v = a[i + j + len / 2] * w;
  233.  
  234.                 a[i + j] = u + v;
  235.                 a[i + j + len / 2] = u - v;
  236.  
  237.                 w *= wlen;
  238.             }
  239.         }
  240.     }
  241.  
  242.     if (invert) {
  243.         for (cd &x : a) x /= n;
  244.     }
  245. }
  246.  
  247. // 2.2) FFT Convolution (polynomial multiplication with real numbers).
  248. //      Multiplies two polynomials using FFT with complex numbers.
  249. //      Parameters:
  250. //        - a: vector of integers (coefficients of the first polynomial)
  251. //        - b: vector of integers (coefficients of the second polynomial)
  252. //      Returns:
  253. //        - A vector of long long integers (rounded from the complex result).
  254. //      Time complexity: O(n log n), where n is the power-of-two size.
  255. //      Important constraints:
  256. //        - The result is rounded to the nearest integer. Precision errors
  257. //          can occur if the coefficients are very large ( > 1e9 ) or the
  258. //          polynomial degree is very high.
  259. //        - Use this when you don't have an NTT-friendly modulus or when
  260. //          you need the exact integer result (not modulo).
  261. vector<long long> convolutionFFT(const vector<int>& a, const vector<int>& b) {
  262.     int n = a.size(), m = b.size();
  263.     if (!n || !m) return {};
  264.  
  265.     int sz = 1;
  266.     while (sz < n + m - 1) sz <<= 1;
  267.  
  268.     vector<cd> fa(a.begin(), a.end()), fb(b.begin(), b.end());
  269.     fa.resize(sz);
  270.     fb.resize(sz);
  271.  
  272.     fft(fa, false);
  273.     fft(fb, false);
  274.  
  275.     for (int i = 0; i < sz; i++) {
  276.         fa[i] *= fb[i];
  277.     }
  278.  
  279.     fft(fa, true);
  280.  
  281.     vector<long long> res(n + m - 1);
  282.     const double EPS = 1e-9;
  283.     for (int i = 0; i < n + m - 1; i++) {
  284.         res[i] = (long long)round(fa[i].real() + EPS);
  285.     }
  286.  
  287.     return res;
  288. }
  289.  
  290. // ===================================================================
  291. // 3) Common Problems Solved with FFT/NTT
  292. //    These are patterns that frequently appear in ECPC/ACPC problems.
  293. //    The idea is to convert the problem into a convolution.
  294. // ===================================================================
  295.  
  296. // 3.1) Counting all pair sums.
  297. //      Given an array of integers, count how many pairs (i, j) have a
  298. //      sum equal to each possible value.
  299. //      Parameters:
  300. //        - arr: vector of integers (the input array)
  301. //      Returns:
  302. //        - A vector 'res' where res[s] = number of ordered pairs with sum = s.
  303. //      Time complexity: O(n log n), where n = max_value - min_value.
  304. //      Important constraints:
  305. //        - The array values must be non-negative. If they can be negative,
  306. //          shift them by the minimum value to make them non-negative.
  307. //        - The result counts ordered pairs (i, j) including i=j.
  308. //          If you need unordered pairs (i < j), adjust the result.
  309. vector<long long> countPairSums(const vector<int>& arr) {
  310.     if (arr.empty()) return {};
  311.  
  312.     int minVal = *min_element(arr.begin(), arr.end());
  313.     int maxVal = *max_element(arr.begin(), arr.end());
  314.  
  315.     int shift = -minVal;
  316.     int size = maxVal - minVal + 1;
  317.  
  318.     vector<int> freq(size, 0);
  319.     for (int x : arr) {
  320.         freq[x + shift]++;
  321.     }
  322.  
  323.     vector<long long> conv = convolutionFFT(freq, freq);
  324.     vector<long long> result(2 * size - 1, 0);
  325.     for (int s = 0; s < (int)conv.size(); s++) {
  326.         int sum = s - 2 * shift;
  327.         if (0 <= sum && sum < (int)result.size()) {
  328.             result[sum] = conv[s];
  329.         }
  330.     }
  331.     return result;
  332. }
  333.  
  334. // 3.2) Counting all pair differences.
  335. //      Given an array of integers, count how many ordered pairs (i, j)
  336. //      have a difference arr[i] - arr[j] equal to each possible value.
  337. //      Parameters:
  338. //        - arr: vector of integers (the input array)
  339. //      Returns:
  340. //        - A vector 'res' where res[d] = number of ordered pairs with difference = d.
  341. //      Time complexity: O(n log n), where n = max_value - min_value.
  342. vector<long long> countPairDifferences(const vector<int>& arr) {
  343.     if (arr.empty()) return {};
  344.  
  345.     int minVal = *min_element(arr.begin(), arr.end());
  346.     int maxVal = *max_element(arr.begin(), arr.end());
  347.  
  348.     int shift = -minVal;
  349.     int size = maxVal - minVal + 1;
  350.  
  351.     vector<int> freq(size, 0);
  352.     for (int x : arr) {
  353.         freq[x + shift]++;
  354.     }
  355.  
  356.     vector<int> revFreq = freq;
  357.     reverse(revFreq.begin(), revFreq.end());
  358.  
  359.     vector<long long> conv = convolutionFFT(freq, revFreq);
  360.     vector<long long> result(2 * size - 1, 0);
  361.     for (int idx = 0; idx < (int)conv.size(); idx++) {
  362.         int diff = idx - (size - 1);
  363.         if (0 <= diff && diff < (int)result.size()) {
  364.             result[diff] = conv[idx];
  365.         }
  366.     }
  367.     return result;
  368. }
  369.  
  370. // 3.3) Counting all subarray sums (for non‑negative arrays).
  371. //      Given an array of non‑negative integers, count how many subarrays
  372. //      have each possible sum.
  373. //      Parameters:
  374. //        - arr: vector of non‑negative integers.
  375. //      Returns:
  376. //        - A vector 'res' where res[s] = number of subarrays with sum = s.
  377. //      Time complexity: O(T log T) where T = total sum of the array.
  378. //      Important constraints:
  379. //        - All elements MUST be non‑negative. The method does not work
  380. //          with negative numbers because prefix sums are not monotonic.
  381. vector<long long> countSubarraySums(const vector<int>& arr) {
  382.     int n = arr.size();
  383.     if (n == 0) return {};
  384.  
  385.     // Check non‑negativity.
  386.     for (int x : arr) {
  387.         if (x < 0) return {}; // Not supported.
  388.     }
  389.  
  390.     int totalSum = accumulate(arr.begin(), arr.end(), 0);
  391.  
  392.     vector<int> prefFreq(totalSum + 1, 0);
  393.     prefFreq[0] = 1; // empty prefix
  394.     int pref = 0;
  395.     for (int x : arr) {
  396.         pref += x;
  397.         prefFreq[pref]++;
  398.     }
  399.  
  400.     // We need sum_{i} freq[i] * freq[i+S] for each S.
  401.     // This is cross‑correlation. Let revFreq[j] = freq[totalSum - j].
  402.     // Then (freq * revFreq)[totalSum - S] = sum_i freq[i] * freq[i+S].
  403.     vector<int> revFreq = prefFreq;
  404.     reverse(revFreq.begin(), revFreq.end());
  405.  
  406.     vector<long long> conv = convolutionFFT(prefFreq, revFreq);
  407.     vector<long long> result(totalSum + 1, 0);
  408.  
  409.     for (int S = 0; S <= totalSum; S++) {
  410.         int idx = totalSum - S;
  411.         if (0 <= idx && idx < (int)conv.size()) {
  412.             result[S] = conv[idx];
  413.         }
  414.         // For S = 0, result[0] = sum_i freq[i]^2 (ordered pairs including i=j).
  415.         // Number of subarrays with sum 0 is sum_i C(freq[i], 2) = (sum_i freq[i]^2 - (n+1)) / 2.
  416.         // We adjust here.
  417.         if (S == 0) {
  418.             long long ordered = result[0];
  419.             long long totalPref = n + 1; // number of prefix sums
  420.             result[0] = (ordered - totalPref) / 2;
  421.         }
  422.     }
  423.  
  424.     return result;
  425. }
  426.  
  427. // 3.4) String matching with wildcards (e.g., '*' matches any character).
  428. //      Given a text string and a pattern string that may contain wildcards,
  429. //      find all positions in the text where the pattern matches.
  430. //      Parameters:
  431. //        - text: the text string (lowercase letters)
  432. //        - pattern: the pattern string (lowercase letters and '*' wildcard)
  433. //      Returns:
  434. //        - A vector of indices (0-based) where the pattern matches.
  435. //      Time complexity: O((n+m) log (n+m)) * alphabet_size
  436. //      Important constraints:
  437. //        - The strings should contain only lowercase letters and '*'.
  438. //        - The wildcard '*' matches any single character.
  439. vector<int> wildcardMatching(const string& text, const string& pattern) {
  440.     int n = text.size(), m = pattern.size();
  441.     if (m > n) return {};
  442.  
  443.     const int ALPHA = 26;
  444.     vector<int> matches(n - m + 1, 0);
  445.  
  446.     for (char c = 'a'; c <= 'z'; c++) {
  447.         vector<int> A(n, 0), B(m, 0);
  448.         for (int i = 0; i < n; i++) {
  449.             if (text[i] == c) A[i] = 1;
  450.         }
  451.         for (int j = 0; j < m; j++) {
  452.             if (pattern[j] == c || pattern[j] == '*') B[j] = 1;
  453.         }
  454.  
  455.         reverse(B.begin(), B.end());
  456.         vector<long long> conv = convolutionFFT(A, B);
  457.  
  458.         for (int i = 0; i <= n - m; i++) {
  459.             matches[i] += conv[i + m - 1];
  460.         }
  461.     }
  462.  
  463.     int required = 0;
  464.     for (char ch : pattern) {
  465.         if (ch != '*') required++;
  466.     }
  467.  
  468.     vector<int> result;
  469.     for (int i = 0; i <= n - m; i++) {
  470.         if (matches[i] == required) result.push_back(i);
  471.     }
  472.     return result;
  473. }
  474.  
  475. // 3.5) Multiplying large integers (Big Integer multiplication).
  476. //      Given two large integers as strings, multiply them.
  477. //      Parameters:
  478. //        - a: string representing the first integer
  479. //        - b: string representing the second integer
  480. //      Returns:
  481. //        - A string representing the product.
  482. //      Time complexity: O(n log n), where n = max(len(a), len(b)).
  483. //      Important constraints:
  484. //        - The input strings should only contain digits (0-9).
  485. string multiplyBigIntegers(const string& a, const string& b) {
  486.     if (a == "0" || b == "0") return "0";
  487.  
  488.     int n = a.size(), m = b.size();
  489.     vector<int> A(n), B(m);
  490.  
  491.     for (int i = 0; i < n; i++) A[i] = a[n - 1 - i] - '0';
  492.     for (int i = 0; i < m; i++) B[i] = b[m - 1 - i] - '0';
  493.  
  494.     vector<long long> conv = convolutionFFT(A, B);
  495.  
  496.     vector<int> res(conv.size() + 1, 0);
  497.     for (int i = 0; i < (int)conv.size(); i++) {
  498.         res[i] += conv[i];
  499.         res[i + 1] += res[i] / 10;
  500.         res[i] %= 10;
  501.     }
  502.  
  503.     while (res.size() > 1 && res.back() == 0) res.pop_back();
  504.  
  505.     string ans;
  506.     for (int i = (int)res.size() - 1; i >= 0; i--) {
  507.         ans.push_back(char('0' + res[i]));
  508.     }
  509.     return ans;
  510. }
  511.  
  512. // ===================================================================
  513. // 4) Advanced Techniques (Formal Power Series)
  514. //    These are more sophisticated operations that appear in harder problems.
  515. // ===================================================================
  516.  
  517. // 4.1) Online NTT (Divide and Conquer DP optimization).
  518. //      Used when dp[i] depends on previous dp values through a convolution.
  519. //      Example: dp[i] = sum_{j < i} dp[j] * f[i-j].
  520. //      This can be computed in O(n log^2 n) using divide and conquer + NTT.
  521. //      Parameters:
  522. //        - dp: vector to be filled (dp[0] should be initialized)
  523. //        - f: the convolution kernel (f[0] is usually 0)
  524. //        - n: number of terms to compute
  525. //      Returns:
  526. //        - Nothing (dp is modified in-place).
  527. //      Time complexity: O(n log^2 n)
  528. //      Important constraints:
  529. //        - dp[0] must be set before calling this function.
  530. //        - The result is computed modulo MOD.
  531. void onlineNTT(vector<int>& dp, const vector<int>& f, int n) {
  532.     function<void(int,int)> solve = [&](int l, int r) {
  533.         if (l == r) return;
  534.         int mid = (l + r) / 2;
  535.         solve(l, mid);
  536.  
  537.         int len1 = mid - l + 1;
  538.         int len2 = r - mid;
  539.         vector<int> A(len1), B(len1 + len2 - 1);
  540.  
  541.         for (int i = l; i <= mid; i++) A[i - l] = dp[i];
  542.         for (int i = 0; i < len1 + len2 - 1; i++) {
  543.             B[i] = (i < (int)f.size() ? f[i] : 0);
  544.         }
  545.  
  546.         vector<int> C = convolution(A, B);
  547.         for (int i = mid + 1; i <= r; i++) {
  548.             dp[i] = (dp[i] + C[i - l]) % MOD;
  549.         }
  550.  
  551.         solve(mid + 1, r);
  552.     };
  553.  
  554.     solve(0, n - 1);
  555. }
  556.  
  557. // 4.2) Polynomial Inverse (formal power series inverse).
  558. //      Given a polynomial A(x), compute its inverse modulo x^n.
  559. //      That is, find B(x) such that A(x) * B(x) ≡ 1 (mod x^n).
  560. //      Parameters:
  561. //        - a: vector of coefficients of A (a[0] must be non‑zero)
  562. //        - n: the number of terms to compute
  563. //      Returns:
  564. //        - A vector of length n representing the inverse polynomial.
  565. //      Time complexity: O(n log n)
  566. //      Important constraints:
  567. //        - a[0] must be invertible modulo MOD.
  568. //        - The result is computed modulo MOD.
  569. vector<int> polynomialInverse(const vector<int>& a, int n) {
  570.     vector<int> res(1, modpow(a[0], MOD - 2, MOD));
  571.     int cur = 1;
  572.  
  573.     while (cur < n) {
  574.         int next = min(cur * 2, n);
  575.         vector<int> f(a.begin(), a.begin() + min((int)a.size(), next));
  576.         vector<int> g = res;
  577.  
  578.         f.resize(next);
  579.         g.resize(next);
  580.  
  581.         vector<int> fg = convolution(f, g);
  582.         fg.resize(next);
  583.  
  584.         for (int i = 0; i < next; i++) {
  585.             fg[i] = (MOD - fg[i]) % MOD;
  586.         }
  587.         fg[0] = (fg[0] + 2) % MOD;
  588.  
  589.         res = convolution(g, fg);
  590.         res.resize(next);
  591.         cur = next;
  592.     }
  593.  
  594.     res.resize(n);
  595.     return res;
  596. }
  597.  
  598. // 4.3) Polynomial Logarithm (log of a formal power series).
  599. //      Computes log(A(x)) modulo x^n.
  600. //      Parameters:
  601. //        - a: vector of coefficients of A (a[0] must be 1)
  602. //        - n: the number of terms to compute
  603. //      Returns:
  604. //        - A vector of length n representing log(A(x)).
  605. //      Time complexity: O(n log n)
  606. //      Important constraints:
  607. //        - a[0] must be 1.
  608. //        - The result is computed modulo MOD.
  609. vector<int> polynomialLog(const vector<int>& a, int n) {
  610.     vector<int> der(max(0, (int)a.size() - 1));
  611.     for (int i = 1; i < (int)a.size(); i++) {
  612.         der[i - 1] = (long long)a[i] * i % MOD;
  613.     }
  614.  
  615.     vector<int> inv = polynomialInverse(a, n);
  616.     vector<int> prod = convolution(der, inv);
  617.     prod.resize(n);
  618.  
  619.     vector<int> res(n, 0);
  620.     for (int i = 1; i < n; i++) {
  621.         res[i] = (long long)prod[i - 1] * modpow(i, MOD - 2, MOD) % MOD;
  622.     }
  623.     return res;
  624. }
  625.  
  626. // 4.4) Polynomial Exponential (exp of a formal power series).
  627. //      Computes exp(A(x)) modulo x^n.
  628. //      Parameters:
  629. //        - a: vector of coefficients of A (a[0] must be 0)
  630. //        - n: the number of terms to compute
  631. //      Returns:
  632. //        - A vector of length n representing exp(A(x)).
  633. //      Time complexity: O(n log n)
  634. //      Important constraints:
  635. //        - a[0] must be 0.
  636. vector<int> polynomialExp(const vector<int>& a, int n) {
  637.     vector<int> res(1, 1);
  638.     int cur = 1;
  639.  
  640.     while (cur < n) {
  641.         int next = min(cur * 2, n);
  642.         vector<int> logRes = polynomialLog(res, next);
  643.         vector<int> diff(next, 0);
  644.  
  645.         for (int i = 0; i < next; i++) {
  646.             diff[i] = (a[i] - logRes[i] + MOD) % MOD;
  647.         }
  648.         diff[0] = (diff[0] + 1) % MOD;
  649.  
  650.         res = convolution(res, diff);
  651.         res.resize(next);
  652.         cur = next;
  653.     }
  654.  
  655.     res.resize(n);
  656.     return res;
  657. }
  658.  
  659. // 4.5) Polynomial Power (raising a polynomial to a power).
  660. //      Computes A(x)^k modulo x^n efficiently.
  661. //      Parameters:
  662. //        - a: vector of coefficients of A(x).
  663. //        - k: the exponent (can be a long long).
  664. //        - n: the number of terms to compute.
  665. //      Returns:
  666. //        - A vector of length n representing A(x)^k.
  667. //      Time complexity: O(n log n)
  668. //      Important constraints:
  669. //        - If a[0] != 0, uses exp(k * log(A)).
  670. //        - If a[0] == 0, it shifts the polynomial to make the first term non-zero.
  671. vector<int> polynomialPower(const vector<int>& a, long long k, int n) {
  672.     int shift = 0;
  673.     while (shift < (int)a.size() && a[shift] == 0) shift++;
  674.  
  675.     if (shift == (int)a.size()) {
  676.         return vector<int>(n, 0);
  677.     }
  678.  
  679.     if ((long long)shift * k >= n) {
  680.         return vector<int>(n, 0);
  681.     }
  682.  
  683.     vector<int> b(a.begin() + shift, a.end());
  684.     int leading = b[0];
  685.  
  686.     int invLeading = modpow(leading, MOD - 2, MOD);
  687.     for (int &x : b) {
  688.         x = (long long)x * invLeading % MOD;
  689.     }
  690.  
  691.     vector<int> logB = polynomialLog(b, n - shift * k);
  692.     for (int &x : logB) {
  693.         x = (long long)x * (k % MOD) % MOD;
  694.     }
  695.  
  696.     vector<int> expLog = polynomialExp(logB, n - shift * k);
  697.  
  698.     int leadingPow = modpow(leading, k, MOD);
  699.     vector<int> res(n, 0);
  700.     for (int i = 0; i < (int)expLog.size(); i++) {
  701.         res[i + shift * k] = (long long)expLog[i] * leadingPow % MOD;
  702.     }
  703.  
  704.     return res;
  705. }
  706.  
  707. // 4.6) Polynomial Square Root.
  708. //      Computes sqrt(A(x)) modulo x^n.
  709. //      Parameters:
  710. //        - a: vector of coefficients of A(x). a[0] must be a quadratic residue.
  711. //        - n: the number of terms to compute.
  712. //      Returns:
  713. //        - A vector of length n representing sqrt(A(x)).
  714. //      Time complexity: O(n log n)
  715. //      Important constraints:
  716. //        - Assumes a[0] = 1 (most common case). For a general a[0], you need
  717. //          to compute its square root modulo MOD.
  718. vector<int> polynomialSqrt(const vector<int>& a, int n) {
  719.     const int INV2 = (MOD + 1) / 2; // modular inverse of 2
  720.     vector<int> res(1, 1); // sqrt(1) = 1
  721.     int cur = 1;
  722.  
  723.     while (cur < n) {
  724.         int next = min(cur * 2, n);
  725.         vector<int> f(a.begin(), a.begin() + min((int)a.size(), next));
  726.         f.resize(next);
  727.  
  728.         vector<int> invRes = polynomialInverse(res, next);
  729.  
  730.         vector<int> prod = convolution(f, invRes);
  731.         prod.resize(next);
  732.  
  733.         for (int i = 0; i < next; i++) {
  734.             res[i] = (long long)(res[i] + prod[i]) * INV2 % MOD;
  735.         }
  736.         res.resize(next);
  737.         cur = next;
  738.     }
  739.  
  740.     res.resize(n);
  741.     return res;
  742. }
  743.  
  744. // ===================================================================
  745. // 5) Utility functions for working with polynomials.
  746. // ===================================================================
  747.  
  748. // 5.1) Trim a polynomial (remove trailing zeros).
  749. vector<int> trimPoly(const vector<int>& a) {
  750.     vector<int> res = a;
  751.     while (!res.empty() && res.back() == 0) res.pop_back();
  752.     if (res.empty()) res.push_back(0);
  753.     return res;
  754. }
  755.  
  756. // 5.2) Evaluate a polynomial at a point x.
  757. long long evalPoly(const vector<int>& a, long long x) {
  758.     long long res = 0;
  759.     for (int i = (int)a.size() - 1; i >= 0; i--) {
  760.         res = (res * x + a[i]) % MOD;
  761.     }
  762.     return res;
  763. }
  764.  
  765. // 5.3) Derivative of a polynomial.
  766. vector<int> derivative(const vector<int>& a) {
  767.     if (a.size() <= 1) return {0};
  768.     vector<int> res(a.size() - 1);
  769.     for (int i = 1; i < (int)a.size(); i++) {
  770.         res[i - 1] = (long long)a[i] * i % MOD;
  771.     }
  772.     return res;
  773. }
  774.  
  775. // 5.4) Integral of a polynomial (constant term is 0).
  776. vector<int> integral(const vector<int>& a) {
  777.     vector<int> res(a.size() + 1, 0);
  778.     for (int i = 0; i < (int)a.size(); i++) {
  779.         res[i + 1] = (long long)a[i] * modpow(i + 1, MOD - 2, MOD) % MOD;
  780.     }
  781.     return res;
  782. }
  783.  
  784. // ===================================================================
  785. // 6) Summary of when to use FFT vs NTT
  786. //    - Use FFT (complex numbers) when:
  787. //        * You need the exact integer result (not modulo).
  788. //        * The coefficients are small enough to avoid precision errors.
  789. //        * The modulus is not NTT-friendly.
  790. //    - Use NTT (modular arithmetic) when:
  791. //        * You need the result modulo a prime.
  792. //        * The modulus is NTT-friendly (e.g., 998244353).
  793. //        * You want exact results without precision issues.
  794. //    - Use Arbitrary Modulus Convolution when:
  795. //        * The modulus is not NTT-friendly but you need exact results.
  796. //        * The modulus is up to ~1e9.
  797. // ===================================================================
  798.  
  799. // ===================================================================
  800. // main() with example usage (you can ignore this part)
  801. // ===================================================================
  802.  
  803. int main() {
  804.     ios::sync_with_stdio(false);
  805.     cin.tie(nullptr);
  806.  
  807.     // Example 1: Polynomial multiplication.
  808.     vector<int> a = {1, 2, 3}; // 1 + 2x + 3x^2
  809.     vector<int> b = {4, 5};    // 4 + 5x
  810.     vector<int> c = convolution(a, b);
  811.     // Expected: 4 + 13x + 22x^2 + 15x^3
  812.  
  813.     cout << "Convolution result: ";
  814.     for (int x : c) cout << x << " ";
  815.     cout << "\n";
  816.  
  817.     // Example 2: Big integer multiplication.
  818.     string bigA = "123456789";
  819.     string bigB = "987654321";
  820.     string product = multiplyBigIntegers(bigA, bigB);
  821.     cout << bigA << " * " << bigB << " = " << product << "\n";
  822.  
  823.     // Example 3: Wildcard matching.
  824.     string text = "abcde";
  825.     string pattern = "a*e";
  826.     vector<int> matches = wildcardMatching(text, pattern);
  827.     cout << "Wildcard matches at positions: ";
  828.     for (int pos : matches) cout << pos << " ";
  829.     cout << "\n";
  830.  
  831.     return 0;
  832. }
Success #stdin #stdout 0s 5316KB
comments (?)
stdin 
Standard input is empty
stdout 
Convolution result: 4 13 22 15 
123456789 * 987654321 = 121932631112635269
Wildcard matches at positions: 0 2
https://ideone.com/1IkGJg
language:
C++ (gcc 8.3)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!