コンパクトに一般modの逆元を求める

ユークリッドの互除法(a,b)->(b,a-floor(a/b)*b) という変換を繰り返すアルゴリズムである。この変換は行列 A={{0,1},{1,-floor(a/b)}}を用いて (a',b')=A(a,b)と書ける。これを繰り返し行うため、互除法の過程は(1,0)=...DCBA(a,b)のような行列積で表せる。

(a,M) で互除法を走らせる。互除法の過程は mod M の上で (1,0)=...DCBA(a,M)=...DCBA(a,0) となる。...DCBAの(1,1)要素を見れば逆元が得られる。

int modinv(int a, int b=mod, int x=1, int y=0) {
	if (b==0) return x<0 ? x+mod : x;
	return modinv(b, a-a/b*b, y, x-a/b*y);
}

重心分解?

2017/08/02(修正) HL分解と重心分解の定義を誤解してた?

HL分解:size[v]<size[c]*2の辺を選択
重心分解:argmax size[c]の辺を選択

いままで重心分解で書いてたけど、HL分解の方が実装が楽?(重心分解はmaxの更新が面倒だと思う)。

struct HLD {
	vector<int> parent, head, vid;

	HLD(const vector<vector<int>> &g) : parent(g.size()), head(g.size()), vid(g.size()) {
		int k = 0;
		vector<int> s(g.size(), 1);
		function<void(int, int)> dfs = [&](int u, int p) {
			for (int v : g[u]) if (v != p) dfs(v, u), s[u] += s[v];
		};
		function<void(int, int, int)> dfs2 = [&](int u, int p, int h) {
			parent[u] = p; head[u] = h; vid[u] = k++;
			for (int v : g[u]) if (v != p && s[u] < s[v] * 2) dfs2(v, u, h);
			for (int v : g[u]) if (v != p && s[u] >= s[v] * 2) dfs2(v, u, v);
		};
		dfs(0, -1);
		dfs2(0, -1, 0);
	}

	template<class T>
	void foreach(int u, int v, T f) {
		for (;; v = parent[v]) {
			if (vid[u] > vid[v]) swap(u, v);
			f(max(vid[u], vid[head[v]]), vid[v]);
			if (head[u] == head[v]) break;
		}
	}

	int lca(int u, int v) {
		for (; head[u] != head[v]; v = parent[head[v]]) if (vid[u] > vid[v]) swap(u, v);
		return vid[u] < vid[v] ? u : v;
	}
};

[1]ではheavy-edgeをsize[v]<size[c]*2と定義している。heavy-light decomposition という表現は見当たらないが、恐らくheavy-light decompositonの由来だろう。

[2]ではcentroid path decomposition on treeとして最大の部分木を繋ぐ辺を heavy-edge としていた。

どちらもheavy-edgeという用語を使っているため、HL分解と一括りにしても問題はない?

参考文献

[1] D. D. Sleator and R. E. Tarjan, A data structure for dynamic trees, in Proc. Thirteenth Annual ACM Symp. on Theory of Computing, 1981.
[2] Paolo Ferragina, Roberto Grossi, Ankur Gupta, Rahul Shah, Jeffrey Scott Vitter, On Searching Compressed String Collections Cache Obliviously, PODS, 2008

コンパクトにmod逆元を求める

mod逆元を求めたいだけなのに繰り返し二乗法を書くのは面倒。そういうときに使える手法。素数mod限定。

long long modinv(long long n) {
	return n == 1 ? n : modinv(MOD % n) * (MOD - MOD / n) % MOD;
}

上の計算式は M%n+⌊M/n⌋n=M≡0 (mod M) に基づいている(除法の原理)。mod 109+7 では再帰の深さが49以下であることを確認している。

計算量解析

雰囲気だけ書いておく。証明にはなっていないし、証明が存在するかどうかも知らない。

M%n<n/2 を満たす n の割合をコンピュータで計算してみると、約 2-log4≈0.6137 だと分かる。2-log4 の確率で値が半分になると仮定すると、再帰の深さの期待値は lg n/(2-log4) になる。実際、mod 109+7での再帰の深さの最大値が49であるのに対して期待値は48.7161なので仮定は正しそうである。よって期待計算量は O(log n) と言えそうである。

yukicoder 550: 夏休みの思い出(1)

https://yukicoder.me/problems/no/550

解法

x3+ax2+bx+c=0 を解く前に x3+ax2+bx+c≡0 (mod 33331) で解く。これは全探索で解くことができ、解の個数が 3 個以下であることが保証されている。

得られた解を α,β,γ とすると、元の方程式の解は α+33331k,β+33331k,γ+33331k の形で表せる。解の範囲が分かっているため、こちらも全探索できる。

#include <iostream>
#include <set>

int main() {
    long long a, b, c;
    std::cin >> a >> b >> c;

    constexpr int MOD = 33331;
    std::set<long long> ans;
    for (__int128 i = 0; i < MOD; i++) {
        if ((i*i*i + a*i*i + b*i + c) % MOD == 0) {
            for (__int128 j = i - MOD*MOD; j <= MOD*MOD; j += MOD) {
                if (j*j*j + a*j*j + b*j + c == 0) ans.insert(j);
            }
        }
    }

    for (long long x : ans) {
        std::cout << x << ' ';
    }
}

python3 (153 bytes)

a,b,c=map(int,input().split())
P=33331
f=lambda x:(x+a)*x*x+x*b+c
print(*sorted(j for i in range(P) if f(i)%P==0 for j in range(i-P*P,P*P,P) if f(j)==0))

解の範囲が狭いことが分かっていれば、より高次の方程式に対しても同じ方法が適用できる。

Codeforces #425 (Div.2) E. Vasya and Shifts

http://codeforces.com/contest/832/problem/E

問題概要

  • 文字列 s1,s2,...,sn がそれぞれ 4 つずつ――全部で4n個――与えられる
  • 使われている文字はabcdeの5種類で、文字列長はすべて m
  • 4n個の文字列の中からいくつか選び、その総和が b になるようなパターン数を求めるクエリが飛んでくる
  • 文字列の和は'a'=0,...,'e'=4として (s+t)[i]=(s[i]+t[i])%5で定義する
  • si を複数回使うとき、使う順序は区別しない。

解法

文字列で考える必要なまったくなくて、単に GF(5) 上の m 次元ベクトルと考えればいい。この問題は

\begin{equation}
\mathbf{b} = c_1 \mathbf{s}_1 + c_2 \mathbf{s}_2 + \cdots +c_n \mathbf{s}_n
\end{equation}

を満たすような \( (c_1,c_2,\ldots,c_n) \) が何通りあるか、という問題に還元できる。これは連立方程式であるから、線形代数の理論を用いると良い。係数行列を A とする。解が存在するならば \(5^{n-\mathrm{rank}A}\) 通り存在する。

したがって、連立方程式が解を持つかどうかを判定できれば良い。\(b\) が与えられるたびに掃き出し法を行うのでは遅いため、\( (A\,\mathbf{b_1}\,\mathbf{b_2}\,\cdots\,\mathbf{b_q}) \)という拡大係数行列を構築しまとめて掃き出し法を行う。

掃き出し法はまあ良いと思うんだけど、解の個数があまり自明ではないので証明のアイディアをつらつらと書いておく。

  • 連立方程式\(Ac=b\) が解 \(Ac'=b\) を持ってるなら、\(A(c-c')=0\) と変形できる
  • 上のように変形できるのであれば、解の個数は連立方程式\(Ax=0\) の解の個数と等しいはず
  • \(Ax=0\)となるような\(x\)全体を行列の\(A\)のカーネルと呼ぶ
  • 解の個数は \( |\mathrm{Ker}A| \) と表せる
  • \( \mathrm{Ker} A \) は線形空間である
  • 次元定理により \( \mathrm{rank}A+\mathrm{dim}\mathrm{Ker} A = \mathrm{dim}V \)――ここで A は V→W の線形写像
  • \(m\) 次元の線形空間は \(m\) 個の基底が存在し、すべての元は \( \lambda_1 u_1 + \lambda_2 u_2 + \cdots + \lambda_m u_m \) の形で一意に定まる
  • \(\mathbb{F}_5\)上の\(m\)次元の線形空間の要素数は \(5^m\) である
#include <iostream>
#include <algorithm>
#include <vector>
#include <string>

constexpr int MOD = 5;

struct modint {
    int n;
    modint(int n = 0) : n(n) {}
};

std::string input() {
    static char buf[1000];
    scanf("%s", buf);
    return std::string(buf);
}

modint inv[5];

modint operator+(modint a, modint b) { return modint((a.n += b.n) >= MOD ? a.n - MOD : a.n); }
modint operator-(modint a, modint b) { return modint((a.n -= b.n) < 0 ? a.n + MOD : a.n); }
modint operator*(modint a, modint b) { return modint(1LL * a.n * b.n % MOD); }
modint operator/(modint a, modint b) { return a * inv[b.n]; }
modint &operator+=(modint &a, modint b) { return a = a + b; }
modint &operator-=(modint &a, modint b) { return a = a - b; }
modint &operator*=(modint &a, modint b) { return a = a * b; }
modint &operator/=(modint &a, modint b) { return a = a / b; }

int main() {
    inv[1] = 1;
    for (int i = 2; i < MOD; i++) {
        inv[i] = inv[MOD % i] * (MOD - MOD / i);
    }

    int n, m;
    std::cin >> n >> m;

    static modint a[500][800];

    for (int i = 0; i < n; i++) {
        std::string s = input();
        for (int j = 0; j < m; j++) {
            a[j][i] = s[j] - 'a';
        }
    }
    int q;
    std::cin >> q;
    for (int i = 0; i < q; i++) {
        std::string s = input();
        for (int j = 0; j < m; j++) {
            a[j][i + n] = s[j] - 'a';
        }
    }

    const int h = m;
    const int w = n + q;
    int r = 0;
    for (int j = 0; j < n && r < h; j++) {
        int p = r;
        for (int i = r; i < h; i++) {
            if (a[i][j].n != 0) {
                p = i;
                break;
            }
        }
        std::swap(a[p], a[r]);
        if (a[r][j].n == 0) continue;
        for (int jj = j + 1; jj < w; jj++) {
            a[r][jj] /= a[r][j];
        }
        a[r][j] = 1;
        for (int ii = r + 1; ii < h; ii++) {
            for (int jj = j + 1; jj < w; jj++) {
                a[ii][jj].n += 25 - a[ii][j].n * a[r][jj].n;
                a[ii][jj].n %= MOD;
            }
            a[ii][j] = 0;
        }
        r++;
    }

    constexpr long long MOD2 = 1e9 + 7;
    long long pw = 1;
    for (int i = 0; i < n - r; i++) {
        pw *= 5;
        pw %= MOD2;
    }

    for (int j = 0; j < q; j++) {
        long long ans = pw;
        for (int i = r; i < m; i++) {
            if (a[i][j + n].n != 0) {
                ans = 0;
            }
        }
        printf("%d\n", (int)ans);
    }
}

187ms。