math

[neal2018]

m = a*(m/a) + m%a

0 = a*(m/a) + m%a Mod m

0 = inv(m%a) * (m/a) + inv(a) Mod m

inv(a) = -inv(m%a) * (m/a) Mod m

inv(a) = inv(m%a) * (m - m/a) Mod m

#include <bits/stdc++.h>
using namespace std;
#define ll long long

ll MOD = 998244353;

ll inv(ll a, ll m) { return a == 1 ? 1 : ((m - m / a) * inv(m % a, m) % m); }

int main(int argc, char const* argv[]) {
    int MAX = 100010;
    vector<ll> f(MAX, 1);
    vector<ll> rf(MAX, 1);
    for (int i = 1; i < MAX; ++i) {
        f[i] = (f[i - 1] * i) % MOD;
        rf[i] = (rf[i - 1] * inv(i, MOD)) % MOD;
    }
    return 0;
}

// alt
vector<ll> f(n + 1, 1), rf(n + 1, 1);
for (int i = 2; i <= n; i++) f[i] = f[i - 1] * i % MOD;
rf[n] = power(f[n], MOD - 2);
for (int i = n - 1; i >= 2; i--) rf[i] = rf[i + 1] * (i + 1) % MOD;

$O(\log(n))$ maybe, open problem

[neal2018]

Fermat's little theorem:

inv(a) = a**(m - 2) Mod m

[neal2018]

mod int class

#include <bits/stdc++.h>
using namespace std;
#define ll long long

constexpr int MOD = 998244353;

constexpr ll norm(ll x)  { return (x % MOD + MOD) % MOD; }
template <class T>
constexpr T power(T a, ll b, T res = 1) {
  for (; b; b /= 2, (a *= a) %= MOD)
    if (b & 1) (res *= a) %= MOD;
  return res;
}
struct Z {
  ll x;
  constexpr Z(ll _x = 0) : x(norm(_x)) {}
  auto operator<=>(const Z &) const = default;
  Z operator-() const { return Z(norm(MOD - x)); }
  Z inv() const { return power(*this, MOD - 2); }
  Z &operator*=(const Z &rhs) { return x = x * rhs.x % MOD, *this; }
  Z &operator+=(const Z &rhs) { return x = norm(x + rhs.x), *this; }
  Z &operator-=(const Z &rhs) { return x = norm(x - rhs.x), *this; }
  Z &operator/=(const Z &rhs) { return *this *= rhs.inv(); }
  Z &operator%=(const ll &rhs) { return x %= rhs, *this; }
  friend Z operator*(Z lhs, const Z &rhs) { return lhs *= rhs; }
  friend Z operator+(Z lhs, const Z &rhs) { return lhs += rhs; }
  friend Z operator-(Z lhs, const Z &rhs) { return lhs -= rhs; }
  friend Z operator/(Z lhs, const Z &rhs) { return lhs /= rhs; }
  friend Z operator%(Z lhs, const ll &rhs) { return lhs %= rhs; }
  friend auto &operator>>(istream &i, Z &z) { return i >> z.x; }
  friend auto &operator<<(ostream &o, const Z &z) { return o << z.x; }
};

modified from https://atcoder.jp/contests/abc220/submissions/26155443

[neal2018]

Euler sieve

vector<int> min_p(maxi + 1), primes;  // minimal prime factor of i; store all primes
for (int i = 2; i <= maxi; i++) {
  if (min_p[i] == 0) {
    min_p[i] = i;
    primes.push_back(i);
  }
  for (auto p : primes) {
    if (i * p > maxi) break;
    min_p[i * p] = p;
    if (i % p == 0) break;
  }
}

line compressed

vector<int> min_p(maxi + 1), primes;
for (int i = 2; i <= maxi; i++) {
  if (min_p[i] == 0) primes.push_back(min_p[i] = i);
  for (auto& p : primes)
    if (i * p > maxi || i % (min_p[i * p] = p) == 0) break;
}

[neal2018]

sqrt decom

for (ll L = 2, R; L <= i; L = R + 1) {
  R = min(i, i / (i / L));
  s = (s + (R - L + 1) * dp[i / L]) % m;
}

[neal2018]

Mobius function

  vector<int> min_p(maxi + 1), mu(maxi + 1), primes;
  mu[1] = 1;
  for (int i = 2; i <= maxi; i++) {
    if (min_p[i] == 0) {
      min_p[i] = i;
      primes.push_back(i);
      mu[i] = -1;
    }
    for (auto p : primes) {
      if (i * p > maxi) break;
      min_p[i * p] = p;
      if (i % p == 0) {
        mu[i * p] = 0;
        break;
      }
      mu[i * p] = -mu[i];
    }
  }

[neal2018]

Möbius function $\mu(n)$

We decompose $$n = \prod_{i=1}^k p_i^{\alpha_i}$$

Then,

$$ \mu(n) = \begin{cases} 0, & \exists \alpha_i \ge 2\\ (-1)^k, & otherwise \end{cases} $$

$$ \varepsilon(n) =\sum_{d|n} \mu(d)=\sum_i \binom{k}{i}(-1)^i = \begin{cases} 1, & n=1\\ 0, & otherwise \end{cases} = [n=1] $$

Dirichlet convolution

$$ (f*g)(n) = \sum{d|n}f(d)g(\frac{n}{d})=\sum{ab=n}f(a)g(b) $$

Identity function

$$ \text{Id}(n)=n $$

Euler's totient function returns the number of integers that are relatively prime and smaller than or equals to $n$

$$ \varphi(n) = \sum_{i\le n} [\text{gcd}(i, n)=1] $$

$$ \sum_{i=1}^n i[\text{gcd}(i, n)=1]= (n\phi(n)+[n=1])/2 $$

Conclusions

$$ \varepsilon * \mu =1 $$

$$ \varphi * 1 =\text{Id} $$

$$ \varphi = \text{Id} *\mu $$

Möbius inversion

$$ f(n)=\sum{d|n}g(d) \Leftrightarrow g(n)=\sum{d|n}\mu(d)f(\frac{n}{d}) $$

$$ f(n)=\sum{n|d}g(d) \Leftrightarrow g(n)=\sum{n|d}\mu(\frac{d}{n})f(d) $$

[neal2018]

Euler's totient function

constexpr int maxi = 5e4 + 10;
vector<int> min_p(maxi + 1), phi(maxi + 1), primes;
phi[1] = 1;
for (int i = 2; i <= maxi; i++) {
  if (min_p[i] == 0) {
    min_p[i] = i;
    primes.push_back(i);
    phi[i] = i - 1;
  }
  for (auto p : primes) {
    if (i * p > maxi) break;
    min_p[i * p] = p;
    if (i % p == 0) {
      phi[i * p] = phi[i] * p;
      break;
    }
    phi[i * p] = phi[i] * phi[p];
  }
}

[neal2018]

Dirichlet prefix sum

a[1]...a[n] -> b[1]...b[n] s.t.

$$ bk = \sum{i|k} a_i $$

for (auto& p : primes)
    for (int i = 1; p * i <= n; i++) cnt[i * p] += cnt[i];

---

Dirichlet suffix sum

a[1]...a[n] -> b[1]...b[n] s.t.

$$ bk = \sum{k|i} a_i $$

for (auto& p : primes)
    for (int i = n/p; i; i--) cnt[i] += cnt[i * p];

[neal2018]

large mod for ntt

constexpr __int128 MOD = 9223372036737335297;

constexpr __int128 norm(__int128 x) { return x < 0 ? (x + MOD) % MOD : x % MOD; }
template <class T>
constexpr T power(T a, __int128 b, T res = 1) {
  for (; b; b /= 2, (a *= a) %= MOD)
    if (b & 1) (res *= a) %= MOD;
  return res;
}
struct Z {
  __int128 x;
  constexpr Z(__int128 _x = 0) : x(norm(_x)) {}
  Z operator-() const { return Z(norm(MOD - x)); }
  Z inv() const { return power(*this, MOD - 2); }
  auto operator<=>(const Z &) const = default;
  Z &operator*=(const Z &rhs) { return x = x * rhs.x % MOD, *this; }
  Z &operator+=(const Z &rhs) { return x = norm(x + rhs.x), *this; }
  Z &operator-=(const Z &rhs) { return x = norm(x - rhs.x), *this; }
  Z &operator/=(const Z &rhs) { return *this *= rhs.inv(); }
  Z &operator%=(const __int128 &rhs) { return x %= rhs, *this; }
  friend Z operator*(Z lhs, const Z &rhs) { return lhs *= rhs; }
  friend Z operator+(Z lhs, const Z &rhs) { return lhs += rhs; }
  friend Z operator-(Z lhs, const Z &rhs) { return lhs -= rhs; }
  friend Z operator/(Z lhs, const Z &rhs) { return lhs /= rhs; }
  friend Z operator%(Z lhs, const __int128 &rhs) { return lhs %= rhs; }
  friend auto operator<=>(Z lhs, const Z &rhs) { return lhs.x <=> rhs.x; }
};

[neal2018]

sqrt decomp

ll n, res = 0;
cin >> n;
for (ll l = 1, r; l <= n; l = r + 1) {
  r = min(n, n / (n / l));
  res += n / l * (r - l + 1);
}
cout << res << "\n";

Mo's algo

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int n, m, maxi = 0;
  cin >> n >> m;
  vector<int> a(n);
  for (auto &x : a) cin >> x, maxi = max(maxi, x);
  vector<array<int, 3>> qs(m);
  for (int i = 0; i < m; i++) {
    auto &[x, y, z] = qs[i];
    cin >> x >> y, x--, z = i;
  }
  int SIZE = max(1, (int)(n / sqrt(m)));
  sort(qs.begin(), qs.end(), [&](auto &x, auto &y) {
    return pair{x[0] / SIZE, x[1]} < pair{y[0] / SIZE, y[1]};
  });
  vector<int> seen(maxi + 1);
  vector<int> res(m);
  int cur = 0;
  auto add = [&](int x) {
    cur -= (seen[x] % 2);
    seen[x]++;
    cur += (seen[x] % 2);
  };
  auto del = [&](int x) {
    cur -= (seen[x] % 2);
    seen[x]--;
    cur += (seen[x] % 2);
  };
  int L = 0, R = 0;  //  store [L, R)
  for (auto &[l, r, id] : qs) {
    while (L > l) add(a[--L]);
    while (R < r) add(a[R++]);
    while (L < l) del(a[L++]);
    while (R > r) del(a[--R]);
    res[id] = cur != 0;
  }
  for (auto &x : res) cout << (x ? "Alice\n" : "Bob\n");
}

Mo's with modification

#include <bits/stdc++.h>
using namespace std;
using ll = long long;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int n, m, maxi = 0, change_cnt = 0, ask_cnt = 0;
  cin >> n >> m;
  vector<int> a(n);
  for (auto &x : a) cin >> x, maxi = max(maxi, x);
  vector<array<int, 4>> qs;
  vector<pair<int, int>> change;
  for (int i = 0, x, y; i < m; i++) {
    string op;
    cin >> op >> x >> y, x--;
    if (op == "R") {
      change.emplace_back(x, y), change_cnt++;
    } else {
      qs.emplace_back(x, y, change_cnt, ask_cnt++);
    }
  }
  int SIZE = (int)cbrt(n * change_cnt) + 1;
  sort(qs.begin(), qs.end(), [&](auto &x, auto &y) {
    return array{x[0] / SIZE, x[1] / SIZE, x[2]} < array{y[0] / SIZE, y[1] / SIZE, y[2]};
  });
  vector<int> seen(maxi + 1);
  vector<int> res(ask_cnt);
  int L = 0, R = 0, now = 0, cur = 0;  //  store [L, R)
  auto add = [&](int x) {
    if (!seen[x]) cur++;
    seen[x]++;
  };
  auto del = [&](int x) {
    seen[x]--;
    if (!seen[x]) cur--;
  };
  auto update = [&](int time) {
    auto [p, v] = change[time];
    if (L <= p && p < R) del(a[p]), add(v);
    change[time].second = a[p], a[p] = v;
  };
  for (auto &[l, r, tt, id] : qs) {
    while (L > l) add(a[--L]);
    while (R < r) add(a[R++]);
    while (L < l) del(a[L++]);
    while (R > r) del(a[--R]);
    while (now < tt) update(now++);
    while (now > tt) update(--now);
    res[id] = cur;
  }
  for (auto &x : res) cout << x << "\n";
}

[neal2018]

int exgcd(int a, int b, int& x, int& y) {
  if (!b) {
    x = 1, y = 0;
    return a;
  }
  int d = exgcd(b, a % b, y, x);
  y -= (a / b) * x;
  return d;
}

[neal2018]

i128 power(i128 a, i128 b, i128 MOD = 1, i128 res = 1) {
  for (; b; b /= 2, (a *= a) %= MOD)
    if (b & 1) (res *= a) %= MOD;
  return res;
}

bool is_prime(ll n) {
  if (n < 2) return false;
  static constexpr int A[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  int s = __builtin_ctzll(n - 1);
  ll d = (n - 1) >> s;
  for (auto a : A) {
    if (a == n) return true;
    ll x = (ll)power(a, d, n);
    if (x == 1 || x == n - 1) continue;
    bool ok = false;
    for (int i = 0; i < s - 1; ++i) {
      x = ll((i128)x * x % n);  // potential overflow!
      if (x == n - 1) {
        ok = true;
        break;
      }
    }
    if (!ok) return false;
  }
  return true;
}

[neal2018]

$$ \beta(k) = \sum_{i=k}^n(-1)^{i-k}\binom{i}{k}\alpha(i) $$

[neal2018]

struct Basis {
  const int SZ = 64;
  vector<ll> a = vector<ll>(SZ);
  Basis(int sz) : SZ(sz), a(sz) {}
  bool insert(ll x) {
    for (int i = SZ - 1; i >= 0; i--) {
      if (x >> i) {
        if (!a[i]) {
          a[i] = x;
          return true;
        } else {
          x ^= a[i];
        }
      }
    }
    return false;
  }
};

[neal2018]

$$ gn = \sum{i=0}^n (-1)^i \binom{n}{i} f_i \Leftrightarrow fn = \sum{i=0}^n(-1)^i\binom{n}{i} g_i $$ $$ gn = \sum{i=0}^n \binom{n}{i} f_i \Leftrightarrow fn = \sum{i=0}^n(-1)^{n-i}\binom{n}{i} g_i $$ $$ gk = \sum{i=k}^n \binom{n}{i} f_i \Leftrightarrow fk = \sum{i=k}^n(-1)^{i-k}\binom{k}{i} g_i $$

[neal2018]

pollard_rho

mt19937_64 rng(chrono::steady_clock::now().time_since_epoch().count());

i128 power(i128 a, i128 b, i128 MOD = 1, i128 res = 1) {
  for (; b; b /= 2, (a *= a) %= MOD)
    if (b & 1) (res *= a) %= MOD;
  return res;
}

bool is_prime(ll n) {
  if (n < 2) return false;
  static constexpr int A[] = {2, 3, 5, 7, 11, 13, 17, 19, 23};
  int s = __builtin_ctzll(n - 1);
  ll d = (n - 1) >> s;
  for (auto a : A) {
    if (a == n) return true;
    ll x = (ll)power(a, d, n);
    if (x == 1 || x == n - 1) continue;
    bool ok = false;
    for (int i = 0; i < s - 1; ++i) {
      x = ll((i128)x * x % n);  // potential overflow!
      if (x == n - 1) {
        ok = true;
        break;
      }
    }
    if (!ok) return false;
  }
  return true;
}

ll pollard_rho(ll x) {
  ll s = 0, t = 0, c = rng() % (x - 1) + 1;
  ll stp = 0, goal = 1, val = 1;
  for (goal = 1;; goal *= 2, s = t, val = 1) {
    for (stp = 1; stp <= goal; ++stp) {
      t = ll(((i128)t * t + c) % x);
      val = ll((i128)val * abs(t - s) % x);
      if ((stp % 127) == 0) {
        ll d = gcd(val, x);
        if (d > 1) return d;
      }
    }
    ll d = gcd(val, x);
    if (d > 1) return d;
  }
}

ll get_max_factor(ll _x) {
  ll max_factor = 0;
  function<void(ll)> fac = [&](ll x) {
    if (x <= max_factor || x < 2) return;
    if (is_prime(x)) {
      max_factor = max_factor > x ? max_factor : x;
      return;
    }
    ll p = x;
    while (p >= x) p = pollard_rho(x);
    while ((x % p) == 0) x /= p;
    fac(x), fac(p);
  };
  fac(_x);
  return max_factor;
}

[neal2018]

matrix inversion

vector<int> operator*(const vector<int>& a, int b) {
  vector<int> res(a.size());
  for (int i = 0; i < a.size(); i++) res[i] = a[i] * b;
  return res;
}

vector<int> operator-(const vector<int>& a, const vector<int>& b) {
  vector<int> res(a.size());
  for (int i = 0; i < a.size(); i++) res[i] = a[i] - b[i];
  return res;
}

vector<int> operator+(const vector<int>& a, int b) {
  vector<int> res(a.size());
  for (int i = 0; i < a.size(); i++) res[i] = a[i] + b;
  return res;
}
vector<int> operator%(const vector<int>& a, int b) {
  vector<int> res(a.size());
  for (int i = 0; i < a.size(); i++) res[i] = a[i] % b;
  return res;
}

vector<vector<int>> inverse(vector<vector<int>> a) {
  int n = int(a.size());
  vector<vector<int>> c(n, vector<int>(n));
  for (int i = 0; i < n; ++i) c[i][i] = 1;
  for (int i = 0; i < n; ++i) {
    for (int j = i; j < n; ++j) {
      if (a[j][i] > 0) {
        swap(a[i], a[j]);
        swap(c[i], c[j]);
        break;
      }
    }
    if (a[i][i] == 0) return {};
    c[i] = c[i] * (1 / a[i][i]);
    a[i] = a[i] * (1 / a[i][i]);
    for (int j = 0; j < n; j++) {
      if (j != i && a[j][i] > 0) {
        c[j] = (c[j] - c[i] * a[j][i] + 2) % 2;
        a[j] = (a[j] - a[i] * a[j][i] + 2) % 2;
      }
    }
  }
  return c;
}

[neal2018]

gaussian elimination

bool is_0(Z v) { return v.x == 0; }
Z abs(Z v) { return v; }
bool is_0(double v) { return abs(v) < 1e-9; }

// 1 => unique solution, 0 => no solution, -1 => multiple solutions
template <typename T>
int gaussian_elimination(vector<vector<T>> &a, int limit) {
    if (a.empty() || a[0].empty()) return -1;
  int h = (int)a.size(), w = (int)a[0].size(), r = 0;
  for (int c = 0; c < limit; c++) {
    int id = -1;
    for (int i = r; i < h; i++) {
      if (!is_0(a[i][c]) && (id == -1 || abs(a[id][c]) < abs(a[i][c]))) {
        id = i;
      }
    }
    if (id == -1) continue;
    if (id > r) {
      swap(a[r], a[id]);
      for (int j = c; j < w; j++) a[id][j] = -a[id][j];
    }
    vector<int> nonzero;
    for (int j = c; j < w; j++) {
      if (!is_0(a[r][j])) nonzero.push_back(j);
    }
    T inv_a = 1 / a[r][c];
    for (int i = r + 1; i < h; i++) {
      if (is_0(a[i][c])) continue;
      T coeff = -a[i][c] * inv_a;
      for (int j : nonzero) a[i][j] += coeff * a[r][j];
    }
    ++r;
  }
  for (int row = h - 1; row >= 0; row--) {
    for (int c = 0; c < limit; c++) {
      if (!is_0(a[row][c])) {
        T inv_a = 1 / a[row][c];
        for (int i = row - 1; i >= 0; i--) {
          if (is_0(a[i][c])) continue;
          T coeff = -a[i][c] * inv_a;
          for (int j = c; j < w; j++) a[i][j] += coeff * a[row][j];
        }
        break;
      }
    }
  } // not-free variables: only it on its line
  for(int i = r; i < h; i++) if(!is_0(a[i][limit])) return 0;   
  return (r == limit) ? 1 : -1;
}

template <typename T>
pair<int,vector<T>> solve_linear(vector<vector<T>> a, const vector<T> &b, int w) {
  int h = (int)a.size();
  for (int i = 0; i < h; i++) a[i].push_back(b[i]);
  int sol = gaussian_elimination(a, w);
  if(!sol) return {0, vector<T>()};
  vector<T> x(w, 0);
  for (int i = 0; i < h; i++) {
    for (int j = 0; j < w; j++) {
      if (!is_0(a[i][j])) {
        x[j] = a[i][w] / a[i][j];
        break;
      }
    }
  }
  return {sol, x};
}