Codeforces 611H New Year and Forgotten Tree

見た目に反してかなり面白い。わからなかったので解説を見たが、行間が異常に広かったので埋めるのに苦労した。この記事ではちゃんと行間を削った解説を書きます。

  • 問題概要

1~Nまでの番号がついた頂点からなる木がある。各辺の両端の頂点番号の桁数のみわかっているとき、木を復元せよ。

  • 解法

頂点を桁数でブロックに分け、各ブロックに1個ずつハブという頂点を設定する。

  • 補題1. 両端が同じ桁数の辺は、すべてそのブロックのハブにつないでよい。
  • 証明. 各ブロック内の辺をどう張ろうと連結成分数は変わらないため。

補題1より、同じ桁数の頂点を結ぶ辺は無視してよい。以下これを無視する。

  • 補題2. 各辺は、少なくとも1個のハブに接続しているとしてよい。
  • 証明. 条件を満たす全域木に、ハブに接続していない辺(i,j)があったとする。iの属するブロックのハブをh_i、 jの属するブロックのハブをh_jとしたとき、辺(i,j)を切断することで、以下の2つの場合のどちらかになる。

場合1. h_iとh_jが異なる連結成分に属する場合、(i,j)を(h_i,h_j)におきかえてよい。

場合2. そうでない場合、(i,j)を(i,h_j)または(h_i,j)のどちらかにおきかえてよい。

  • 補題3. 頂点をハブに限ったグラフは、連結であるとしてよい。すなわち、ハブ同士を結ぶ辺が(ハブ数-1)本あるとしてよい。
  • 証明. そうでないとすると、補題2より、複数のハブに接続するハブでない頂点がある。この頂点をxとし、xのブロックのハブをh_xとする。xを取り除くとグラフはいくつかの連結成分に分かれるが、その中にはh_xを含む連結成分がある。それ以外の連結成分からxにのびる辺を、h_xに付け替えればよい。
  • 補題4. ハブでない頂点からは、ちょうど一つのハブへ辺がはられており、それ以外の辺ははられていない。
  • 証明. 補題3より、ハブ同士を結ぶ辺が(ハブ数-1)本ある。これ以外の辺の本数はハブでない頂点数と等しく、よって連結性と補題2より示される。


以上より、以下のアルゴリズムが得られる。

  • まず、各ブロックでハブを決め、同じブロック内の辺を適切に削除する
  • その後、ハブからなる全域木を全探索する。(個数は6^4=1296個以下)
  • 全域木を固定すると、各々のブロックの2つ組に対してあと何本辺をはる必要があるかが求まり、さらに各ブロックからあと何本の辺が出るかもわかる。補題4より、辺を新しく1本はるごとにハブでない頂点が1個消費されるので、このような辺のはり方が存在するかどうかは、
  • source -> 「ブロックの2つ組に対応する頂点」に容量がそのブロック間にはる必要のある辺の本数

「ブロックの2つ組に対応する頂点」 -> 「ブロックに対応する頂点」に容量無限

「ブロックに対応する頂点」 -> sinkに容量がブロックから出る辺の本数

の辺をはった最大流が必要な辺の本数と等しいかどうかで判定できる。さらに、これが可能な場合、同じく最大流を出力すれば辺のはり方が得られる。



ぼくの実装だとブロックが1個の場合のコーナーケースに引っかかったがそれ以外はバグなくうまくいった。こういうやりたくない実装をそこそこの精度で合わせられるのは成長な気がする。

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<string>
#include<iostream>
using namespace std;
#define SIZE 100
vector<int>pat[SIZE];
vector<int>cap[SIZE];
vector<int>rev[SIZE];
//vector<bool>ise[SIZE];
//vector<bool>isvl[SIZE];
void adde(int a, int b, int c)
{
	pat[a].push_back(b);
	pat[b].push_back(a);
	cap[a].push_back(c);
	cap[b].push_back(0);
	rev[a].push_back(pat[b].size() - 1);
	rev[b].push_back(pat[a].size() - 1);
}
int frv[SIZE], fre[SIZE];
bool flag[SIZE];
void dfs(int node)
{
	if (flag[node])return;
	flag[node] = true;
	for (int i = 0; i<pat[node].size(); i++)
	{
		if (cap[node][i]>0 && (!flag[pat[node][i]]))
		{
			frv[pat[node][i]] = node;
			fre[pat[node][i]] = i;
			dfs(pat[node][i]);
		}
	}
}
int nod;
//int flow=0;
int maxflow(int st, int go)
{
	int ret = 0;
	for (;;)
	{
		fill(flag, flag + nod, false);
		dfs(st);
		if (!flag[go])return ret;
		int mini = 1000000000;
		int now = go;
		for (;;)
		{
			if (now == st)break;
			int t = fre[now];
			now = frv[now];
			mini = min(mini, cap[now][t]);
		}
		ret += mini;
		//flow++;
		now = go;
		for (;;)
		{
			if (now == st)break;
			int t = fre[now];
			int nx = now;
			now = frv[now];
			cap[now][t] -= mini;
			cap[nx][rev[now][t]] += mini;
		}
	}
}
int dat[6][6];
int rem[6];
int pt[6];
int base[6];
typedef pair<int, int>pii;
#define SIZE 6
class unionfind
{
public:
	int par[SIZE];
	int ran[SIZE];
	int ren[SIZE];
	void init()
	{
		for (int i = 0; i<SIZE; i++)
		{
			par[i] = i;
			ran[i] = 0;
			ren[i] = 1;
		}
	}
	int find(int a)
	{
		if (a == par[a])return a;
		else return par[a] = find(par[a]);
	}
	void unite(int a, int b)
	{
		a = find(a);
		b = find(b);
		if (a == b)return;
		if (ran[a]>ran[b])
		{
			par[b] = a;
			ren[a] += ren[b];
		}
		else
		{
			par[a] = b;
			ren[b] += ren[a];
		}
		if (ran[a] == ran[b])ran[b]++;
	}
};
int mxdig;
bool isok = false;
vector<pii>ans;
void calc(vector<pii>vec)
{
	//for (int i = 0; i < vec.size(); i++)printf("%d %d  ", vec[i].first, vec[i].second); printf("\n");
	nod = 100;
	for (int i = 0; i < 100; i++)
	{
		pat[i].clear();
		cap[i].clear();
		rev[i].clear();
	}
	int d[6][6];
	for (int i = 0; i <= mxdig; i++)
	{
		for (int j = 0; j <= mxdig; j++)
		{
			if (i != j)d[i][j] = dat[i][j];
		}
	}
	for (int i = 0; i < vec.size(); i++)d[vec[i].first][vec[i].second]--, d[vec[i].second][vec[i].first]--;
	for (int i = 0; i <= mxdig; i++)
	{
		for (int j = i + 1; j <= mxdig; j++)
		{
			//printf("%d %d  %d\n", i, j, d[i][j]);
			if (d[i][j] < 0)return;
			adde(0, i * 6 + j + 1, d[i][j]);
			adde(i * 6 + j + 1, 50 + i, 1000000000);
			adde(i * 6 + j + 1, 50 + j, 1000000000);
		}
	}
	int mok = 0;
	for (int i = 0; i <= mxdig; i++)
	{
		adde(50 + i, 99, rem[i] - 1);
		mok += rem[i] - 1;
	}
	int fl = maxflow(0, 99);
	if (fl == mok)
	{
		isok = true;
		for (int i = 0; i < vec.size(); i++)ans.push_back(make_pair(base[vec[i].first], base[vec[i].second]));
		for (int i = 0; i <= mxdig; i++)
		{
			for (int j = i + 1; j <= mxdig; j++)
			{
				for (int k = 0; k < pat[i * 6 + j + 1].size(); k++)
				{
					if (pat[i * 6 + j + 1][k] == 50 + i)
					{
						for (int l = 0; l < 1000000000 - cap[i * 6 + j + 1][k]; l++)
						{
							ans.push_back(make_pair(pt[i]++, base[j]));
						}
					}
					if (pat[i * 6 + j + 1][k] == 50 + j)
					{
						for (int l = 0; l < 1000000000 - cap[i * 6 + j + 1][k]; l++)
						{
							ans.push_back(make_pair(base[i], pt[j]++));
						}
					}
				}
			}
		}
	}
}
void dfs(vector<pii>vec, int pt1, int pt2)
{
	if (isok)return;
	unionfind uf;
	uf.init();
	for (int i = 0; i < vec.size(); i++)uf.unite(vec[i].first, vec[i].second);
	int n1 = pt1, n2 = pt2 + 1;
	if (n2 == mxdig + 1)
	{
		n1++;
		n2 = n1 + 1;
	}
	if (n2 == mxdig + 1)
	{
		if (vec.size() == mxdig)
		{
			calc(vec);
		}
		return;
	}
	dfs(vec, n1, n2);
	if (uf.find(n1) != uf.find(n2))
	{
		vec.push_back(make_pair(n1, n2));
		dfs(vec, n1, n2);
	}
}
int main()
{
	int num;
	scanf("%d", &num);
	for (int i = 0; i < num - 1; i++)
	{
		string sa, sb;
		cin >> sa >> sb;
		dat[sa.size() - 1][sb.size() - 1]++;
		if (sa.size() != sb.size())dat[sb.size() - 1][sa.size() - 1]++;
	}
	for (int i = 1; i <= num; i++)
	{
		if (i <= 9)rem[0]++;
		else if (i <= 99)rem[1]++;
		else if (i <= 999)rem[2]++;
		else if (i <= 9999)rem[3]++;
		else if (i <= 99999)rem[4]++;
		else rem[5]++;
	}
	for (int i = 5; i >= 0; i--)
	{
		if (rem[i] > 0)
		{
			mxdig = i;
			break;
		}
	}
	if(mxdig==0)
	{
		for (int i = 0; i < num - 1; i++)printf("1 %d\n", i + 2);
		return 0;
	}
	pt[0] = 1, pt[1] = 10, pt[2] = 100, pt[3] = 1000, pt[4] = 10000, pt[5] = 100000;
	for (int i = 0; i < 6; i++)base[i] = pt[i];
	for (int i = 0; i < 6; i++)
	{
		pt[i]++;
		for (int j = 0; j < dat[i][i]; j++)
		{
			ans.push_back(make_pair(base[i], pt[i]++));
			rem[i]--;
			if (rem[i] <= 0)
			{
				printf("-1\n");
				return 0;
			}
		}
	}
	//for (int i = 0; i < ans.size(); i++)printf("%d %d\n", ans[i].first, ans[i].second);
	vector<pii>zv;
	dfs(zv, 0, 0);
	if (!isok)
	{
		printf("-1\n");
	}
	else
	{
		for (int i = 0; i < ans.size(); i++)printf("%d %d\n", ans[i].first, ans[i].second);
	}
}

Codeforces 633H Fibonacci-ish II

問題概要: 区間をソートした後重複を取り除き、さらにそのi項目にフィボナッチ数列のi項目をかけたものの和を求めるクエリを30000個くらい処理せよ。

解法:

最初に思いついたのが、区間を平方分割して、クエリを終点でソートし、sqrt(N)個のsegtreeを始点をずらして作って一番近いところを少し動かしてクエリにこたえるという方法で、実装したらTLEしたので解説を見た。Moというものがあるらしい。(オーダーは同じ)

Moをやってもなおバケットサイズを気を付けないと通らないしpretest形式で出すものではないと思う。フィボナッチをのせるsegtreeは大昔に一度書いたけど完全に忘却していた。遅延更新が異常な量あるしそりゃ定数倍も重いわという感じ。

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<set>
using namespace std;
#define SIZE 32768
int mod;
int fib[100000];
int rfib[100000];
class BIT
{
public:
	int bit[SIZE + 1];
	void add(int a, int b)
	{
		a++;
		for (;;)
		{
			bit[a] += b;
			a += a&-a;
			if (a > SIZE)return;
		}
	}
	int get(int a)
	{
		a++;
		int ret = 0;
		for (;;)
		{
			ret += bit[a];
			a -= a&-a;
			if (a == 0)return ret;
		}
	}
};
class segtree
{
public:
	int seg[SIZE * 2];
	int pl[SIZE * 2];
	int flag[SIZE * 2];
	void update(int beg, int end, int node, int lb, int ub, int n)
	{
		if (ub < beg || end < lb)return;
		if (beg <= lb&&ub <= end)
		{
			if (n == 1)
			{
				int t = pl[node];
				pl[node] = seg[node];
				seg[node] = (seg[node] + t) % mod;
				flag[node]++;
			}
			else
			{
				int t = seg[node];
				seg[node] = pl[node];
				pl[node] = (t + mod - pl[node]) % mod;
				flag[node]--;
			}
			return;
		}
		if (flag[node])
		{
			if (flag[node] > 0)
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2 + 1] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
			}
			else
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2 + 1] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
			}
			flag[node * 2] += flag[node];
			flag[node * 2 + 1] += flag[node];
			flag[node] = 0;
		}
		update(beg, end, node * 2, lb, (lb + ub) / 2, n);
		update(beg, end, node * 2 + 1, (lb + ub) / 2 + 1, ub, n);
		seg[node] = (seg[node * 2] + seg[node * 2 + 1])%mod;
		pl[node] = (pl[node * 2] + pl[node * 2 + 1])%mod;
	}
	void add(int p, int node, int lb, int ub, int t1, int t2)
	{
		if (ub < p || p < lb)return;
		if (lb == ub)
		{
			seg[node] += t1;
			pl[node] += t2;
			seg[node] %= mod;
			pl[node] %= mod;
		//	printf("add %d %d %d\n", lb, t1, t2);
			return;
		}
		if (flag[node])
		{
			if (flag[node] > 0)
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2 + 1] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
			}
			else
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2 + 1] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
			}
			flag[node * 2] += flag[node];
			flag[node * 2 + 1] += flag[node];
			flag[node] = 0;
		}
		add(p, node * 2, lb, (lb + ub) / 2, t1, t2);
		add(p, node * 2 + 1, (lb + ub) / 2 + 1, ub, t1, t2);
		seg[node] = (seg[node * 2] + seg[node * 2 + 1]) % mod;
		pl[node] = (pl[node * 2] + pl[node * 2 + 1]) % mod;
		//printf("node: %d %d %d\n", node, seg[node], pl[node]);
	}
};
#define B 300
segtree tree;
BIT bi;
typedef pair<int, int>pii;
typedef pair<pii, pii>pi4;
int ans[30000];
int cnt[30000];
int main()
{
	int num;
	scanf("%d%d", &num, &mod);
	vector<int>vec;
	vector<int>zat;
	for (int i = 0; i < num; i++)
	{
		int z;
		scanf("%d", &z);
		vec.push_back(z);
		zat.push_back(z);
	}
	sort(zat.begin(), zat.end());
	int query;
	scanf("%d", &query);
	vector<pi4>dat;
	for (int i = 0; i < query; i++)
	{
		int za, zb;
		scanf("%d%d", &za, &zb);
		za--;
		zb--;
		dat.push_back(make_pair(make_pair(zb / B, za), make_pair(zb, i)));
	}
	fib[1] = rfib[1] = 1;
	for (int i = 2; i < 50000; i++)
	{
		fib[i] = (fib[i - 1] + fib[i - 2]) % mod;
		rfib[i] = (rfib[i - 2] + mod - rfib[i - 1]) % mod;
	}
	sort(dat.begin(), dat.end());
	int nl = 0, nr = -1;
	for (int i = 0; i < query; i++)
	{
		int beg = dat[i].first.second, end = dat[i].second.first, idx = dat[i].second.second;
		for (int p=0;;p++)
		{
			//printf("%d %d %d\n", nl, nr, tree.seg[1]);
			if (nl > beg)
			{
				nl--;
				int low = lower_bound(zat.begin(), zat.end(), vec[nl]) - zat.begin();
				if (cnt[low] == 0)
				{
					int t = bi.get(low);
					tree.add(low, 1, 0, SIZE - 1, vec[nl] % mod*fib[t + 1] % mod, vec[nl] % mod*fib[t] % mod);
					tree.update(low + 1, SIZE - 1, 1, 0, SIZE - 1, 1);
					bi.add(low, 1);
				}
				cnt[low]++;
			}
			else if (nr < end)
			{
				nr++;
				int low = lower_bound(zat.begin(), zat.end(), vec[nr]) - zat.begin();
				if (cnt[low] == 0)
				{
					int t = bi.get(low);
					tree.add(low, 1, 0, SIZE - 1, vec[nr] % mod*fib[t + 1] % mod, vec[nr] % mod*fib[t] % mod);
					tree.update(low + 1, SIZE - 1, 1, 0, SIZE - 1, 1);
					bi.add(low, 1);
				}
				cnt[low]++;
			}
			else if (nl < beg)
			{
				int low = lower_bound(zat.begin(), zat.end(), vec[nl]) - zat.begin();
				if (cnt[low] == 1)
				{
					int t = bi.get(low - 1);
					tree.add(low, 1, 0, SIZE - 1, (mod - vec[nl] % mod*fib[t + 1] % mod) % mod, (mod - vec[nl] % mod*fib[t] % mod) % mod);
					tree.update(low + 1, SIZE - 1, 1, 0, SIZE - 1, -1);
					bi.add(low, -1);
				}
				nl++;
				cnt[low]--;
			}
			else if (nr > end)
			{
				int low = lower_bound(zat.begin(), zat.end(), vec[nr]) - zat.begin();
				if (cnt[low] == 1)
				{
					int t = bi.get(low - 1);
					tree.add(low, 1, 0, SIZE - 1, (mod - vec[nr] % mod*fib[t + 1] % mod) % mod, (mod - vec[nr] % mod*fib[t] % mod) % mod);
					tree.update(low + 1, SIZE - 1, 1, 0, SIZE - 1, -1);
					bi.add(low, -1);
				}
				nr--;
				cnt[low]--;
			}
			else break;
		}
		ans[idx] = tree.seg[1];
	}
	for (int i = 0; i < query; i++)printf("%d\n", ans[i]);
}

以下TLE解

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<set>
using namespace std;
#define SIZE 32768
int mod;
int fib[100000];
int rfib[100000];
class BIT
{
public:
	int bit[SIZE + 1];
	void add(int a, int b)
	{
		a++;
		for (;;)
		{
			bit[a] += b;
			a += a&-a;
			if (a > SIZE)return;
		}
	}
	int get(int a)
	{
		a++;
		int ret = 0;
		for (;;)
		{
			ret += bit[a];
			a -= a&-a;
			if (a == 0)return ret;
		}
	}
};
class segtree
{
public:
	int seg[SIZE * 2];
	int pl[SIZE * 2];
	int flag[SIZE * 2];
	void update(int beg, int end, int node, int lb, int ub, int n)
	{
		if (ub < beg || end < lb)return;
		if (beg <= lb&&ub <= end)
		{
			if (n == 1)
			{
				int t = pl[node];
				pl[node] = seg[node];
				seg[node] = (seg[node] + t) % mod;
				flag[node]++;
			}
			else
			{
				int t = seg[node];
				seg[node] = pl[node];
				pl[node] = (t + mod - pl[node]) % mod;
				flag[node]--;
			}
			return;
		}
		if (flag[node])
		{
			if (flag[node] > 0)
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2 + 1] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
			}
			else
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2 + 1] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
			}
			flag[node * 2] += flag[node];
			flag[node * 2 + 1] += flag[node];
			flag[node] = 0;
		}
		update(beg, end, node * 2, lb, (lb + ub) / 2, n);
		update(beg, end, node * 2 + 1, (lb + ub) / 2 + 1, ub, n);
		seg[node] = (seg[node * 2] + seg[node * 2 + 1])%mod;
		pl[node] = (pl[node * 2] + pl[node * 2 + 1])%mod;
	}
	void add(int p, int node, int lb, int ub, int t1, int t2)
	{
		if (ub < p || p < lb)return;
		if (lb == ub)
		{
			seg[node] += t1;
			pl[node] += t2;
			seg[node] %= mod;
			pl[node] %= mod;
		//	printf("add %d %d %d\n", lb, t1, t2);
			return;
		}
		if (flag[node])
		{
			if (flag[node] > 0)
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*fib[flag[node] + 1] + b*fib[flag[node]]) % mod;
				pl[node * 2 + 1] = (a*fib[flag[node]] + b*fib[flag[node] - 1]) % mod;
			}
			else
			{
				int a = seg[node * 2], b = pl[node * 2];
				seg[node * 2] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
				a = seg[node * 2 + 1], b = pl[node * 2 + 1];
				seg[node * 2 + 1] = (a*rfib[-flag[node] - 1] + b*rfib[-flag[node]]) % mod;
				pl[node * 2 + 1] = (a*rfib[-flag[node]] + b*rfib[-flag[node] + 1]) % mod;
			}
			flag[node * 2] += flag[node];
			flag[node * 2 + 1] += flag[node];
			flag[node] = 0;
		}
		add(p, node * 2, lb, (lb + ub) / 2, t1, t2);
		add(p, node * 2 + 1, (lb + ub) / 2 + 1, ub, t1, t2);
		seg[node] = (seg[node * 2] + seg[node * 2 + 1]) % mod;
		pl[node] = (pl[node * 2] + pl[node * 2 + 1]) % mod;
		//printf("node: %d %d %d\n", node, seg[node], pl[node]);
	}
};
#define B 200
#define NB 150
segtree tree[NB];
BIT bi[NB];
typedef pair<int, int>pii;
vector<pii>dat[30000];
int ans[30000];
int main()
{
	int num;
	scanf("%d%d", &num, &mod);
	vector<int>vec;
	vector<int>zat;
	for (int i = 0; i < num; i++)
	{
		int z;
		scanf("%d", &z);
		vec.push_back(z);
		zat.push_back(z);
	}
	sort(zat.begin(), zat.end());
	int query;
	scanf("%d", &query);
	for (int i = 0; i < query; i++)
	{
		int za, zb;
		scanf("%d%d", &za, &zb);
		za--;
		zb--;
		dat[zb].push_back(make_pair(za, i));
	}
	fib[1] = rfib[1] = 1;
	for (int i = 2; i < 50000; i++)
	{
		fib[i] = (fib[i - 1] + fib[i - 2]) % mod;
		rfib[i] = (rfib[i - 2] + mod - rfib[i - 1]) % mod;
	}
	//for (int i = 0; i < 10; i++)printf(" %d\n", rfib[i]);
	for (int i = 0; i < num; i++)
	{
		//printf("now: %d\n", i);
		for (int j = 0; j < NB; j++)
		{
			if (i >= j*B)
			{
				int low = lower_bound(zat.begin(), zat.end(), vec[i]) - zat.begin();
				if (bi[j].get(low) == bi[j].get(low - 1))
				{
					int t = bi[j].get(low);
					tree[j].add(low, 1, 0, SIZE - 1, vec[i] % mod*fib[t + 1] % mod, vec[i] % mod*fib[t] % mod);
					tree[j].update(low + 1, SIZE - 1, 1, 0, SIZE - 1, 1);
					//printf(" %d  %d %d  %d %d\n",j, low, vec[i], vec[i] % mod*fib[t + 1] % mod, vec[i] % mod*fib[t] % mod);
					bi[j].add(low, 1);
				}
			}
		}
		/*for (int j = 0; j < 5; j++)
		{
			printf("%d ", tree[j].seg[1]);
		}
		printf("\n");*/
		for (int j = 0; j < dat[i].size(); j++)
		{
			int b = (dat[i][j].first + B - 1) / B;
			if (i <= b * B)
			{
				vector<int>v;
				for (int k = dat[i][j].first; k <= i; k++)
				{
					v.push_back(vec[k]);
				}
				sort(v.begin(), v.end());
				int now = -1, pt = 0;
				for (int k = 0; k < v.size(); k++)
				{
					if (now != v[k])
					{
						now = v[k];
						ans[dat[i][j].second] += now%mod*fib[pt + 1];
						pt++;
						ans[dat[i][j].second] %= mod;
					}
				}
				continue;
			}
			//printf("query: %d %d  %d\n", dat[i][j].first, i, b);
			vector<pii>zv;
			for (int k = b*B - 1; k >= dat[i][j].first; k--)
			{
				int low = lower_bound(zat.begin(), zat.end(), vec[k]) - zat.begin();
				if (bi[b].get(low) == bi[b].get(low - 1))
				{
					int t = bi[b].get(low);
					tree[b].add(low, 1, 0, SIZE - 1, vec[k] % mod*fib[t + 1] % mod, vec[k] % mod*fib[t] % mod);
					tree[b].update(low + 1, SIZE - 1, 1, 0, SIZE - 1, 1);
					//printf("%d %d %d, %d\n", low, vec[k] % mod*fib[t + 1] % mod, vec[k] % mod*fib[t] % mod, low + 1);
					zv.push_back(make_pair(vec[k], t));
					bi[b].add(low, 1);
				}
			}
			ans[dat[i][j].second] = tree[b].seg[1];
			for (int k = zv.size() - 1; k >= 0; k--)
			{
				int low = lower_bound(zat.begin(), zat.end(), zv[k].first) - zat.begin();
				tree[b].update(low + 1, SIZE - 1, 1, 0, SIZE - 1, -1);
				tree[b].add(low, 1, 0, SIZE - 1, (mod - zv[k].first % mod*fib[zv[k].second + 1] % mod) % mod, (mod - zv[k].first % mod*fib[zv[k].second] % mod) % mod);
				//printf("%d %d %d, %d\n", low, 0, 0, low + 1);
				bi[b].add(low, -1);
			}
		}
		/*for (int j = 0; j < 5; j++)
		{
			printf("%d ", tree[j].seg[1]);
		}
		printf("\n");*/
	}
	for (int i = 0; i < query; i++)printf("%d\n", ans[i]);
}

Codeforces 668F Little Artem and Graph

http://codeforces.com/contest/668/problem/F

問題概要: サイズkのクリークから始め、サイズkのクリークの全点と新たに結ぶことで頂点を増やしていって作ったグラフの全域木の個数を求めよ。

解法: 解説は木分解とかしてるしantaさんも意味不明なライブラリをはっていて怖いが、普通に行列木するとできる。
最初は非0要素が高々kn個くらいの疎行列で、cliqueに辺をはるので掃き出しても非0要素は増えないことを使って自然に掃き出していくだけ。

コンテスト中に思いついてたのに対角成分に0が生えると逆元が取れなくて詰むとかいう勘違いをしていた。対角成分に0が生えたら答えは0なんじゃ。(つらい)

これ通せてれば2位だったから惜しいなぁ...

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<map>
using namespace std;
typedef long long ll;
ll mod = 1000000007;
ll po(ll a, ll b)
{
	if (b == 0)return 1;
	ll z = po(a, b / 2);
	z = z*z%mod;
	if (b % 2 == 1)z = z*a%mod;
	return z;
}
ll inv(ll a)
{
	return po(a, mod - 2);
}
map<int, ll>ma[10000];
int deg[10000];
vector<int>frm[10000];
int main()
{
	int num, gen;
	scanf("%d%d", &num, &gen);
	for (int i = 0; i < gen; i++)
	{
		for (int j = 0; j < gen; j++)
		{
			if (i != j)ma[i][j] = mod - 1;
			if (i < j)frm[j].push_back(i);
		}
		deg[i] = gen - 1;
	}
	for (int i = gen; i < num; i++)
	{
		for (int j = 0; j < gen; j++)
		{
			int z;
			scanf("%d", &z);
			z--;
			ma[i][z] = ma[z][i] = mod - 1;
			deg[i]++;
			deg[z]++;
			frm[i].push_back(z);
		}
	}
	for (int i = 0; i < num; i++)ma[i][i] = deg[i];
	ll ans = 1;
	for (int i = num - 1; i >= 1; i--)
	{
		ll t = ma[i][i];
		ll iv = inv(mod - t);
		for (int j = 0; j < frm[i].size(); j++)
		{
			ll tim = iv*ma[frm[i][j]][i] % mod;
			for (int k = 0; k < frm[i].size(); k++)
			{
				ma[frm[i][j]][frm[i][k]] = (ma[frm[i][j]][frm[i][k]] + ma[i][frm[i][k]] * tim) % mod;
			}
		}
		ans *= t;
		ans %= mod;
	}
	printf("%I64d\n", ans);
}

Codeforces 674G Choosing Ads

VK Cup Round 3のボス問。

  • 問題概要

区間を埋めるクエリと、区間のうちp%以上を占める数をすべて出力せよというクエリに答えよ。クエリ数, 数の範囲は150000, 20<=p

  • 解法

区間のうち1/5以上を占めるものしか出力しなくていいので、区間からランダムに何回か選べば答えの候補は全部列挙できるはず。
よって、ある位置の値を取得するのと区間を埋めるのができれば答えの候補は全部列挙できる。(ただし、前者はたくさんやるのでlogがつかないように平方分割でやる。)

列挙したら、その区間のその数の個数を数える。これも平方分割でできる。このとき、得られたものすべてを見ていると150000のsqrtにさらに何かがかかって大変なので特定個数以上現れたものを出力するなり上から特定番目のものまでを出力するなりすることにする。

これを素直に(ランダムに選ぶ点の数とか実際に個数を数える基準とかのパラメータを適切に調整して)実装するとTLEするので、個数を数えるクエリの定数倍高速化を頑張る。例えば区間の端のブロックを左か右の小さいほうを見るとかをするとバケットサイズも相まって定数倍1/sqrt(2)とかをつけられたりする。

バケットサイズ等のパラメータをいじり、定数倍高速化を頑張ってやっとAC。パラメータをちゃんと考えれば間違える確率は0に収束するから乱択と言っても一応解法の正当性は証明できるけど、これ絶対想定じゃないよね...

Jun/09/2016 16:47 Compilation error [main tests] → 18342397
Jun/09/2016 16:48 Compilation error [main tests] → 18342410
Jun/09/2016 16:49 Compilation error [main tests] → 18342426
Jun/09/2016 16:51 Compilation error [main tests] → 18342466
Jun/09/2016 16:53 Wrong answer on test 24 [main tests] → 18342487
Jun/09/2016 16:54 Wrong answer on test 31 [main tests] → 18342497
Jun/09/2016 16:55 Wrong answer on test 31 [main tests] → 18342514
Jun/09/2016 16:56 Wrong answer on test 31 [main tests] → 18342531
Jun/09/2016 17:00 Wrong answer on test 31 [main tests] → 18342581
Jun/09/2016 17:02 Wrong answer on test 31 [main tests] → 18342614
Jun/09/2016 17:04 Wrong answer on test 31 [main tests] → 18342651
Jun/09/2016 17:08 Wrong answer on test 31 [main tests] → 18342741
Jun/09/2016 18:00 Time limit exceeded on test 50 [main tests] → 18343573
Jun/09/2016 18:01 Time limit exceeded on test 50 [main tests] → 18343586
Jun/09/2016 18:02 Time limit exceeded on test 50 [main tests] → 18343597
Jun/09/2016 18:05 Time limit exceeded on test 50 [main tests] → 18343627
Jun/09/2016 18:06 Time limit exceeded on test 50 [main tests] → 18343649
Jun/09/2016 18:08 Wrong answer on test 43 [main tests] → 18343668
Jun/09/2016 18:09 Time limit exceeded on test 50 [main tests] → 18343681
Jun/09/2016 18:11 Time limit exceeded on test 50 [main tests] → 18343713
Jun/09/2016 18:13 Wrong answer on test 43 [main tests] → 18343741
Jun/09/2016 18:13 Time limit exceeded on test 50 [main tests] → 18343750
Jun/09/2016 18:17 Time limit exceeded on test 48 [main tests] → 18343806
Jun/09/2016 18:18 Time limit exceeded on test 47 [main tests] → 18343818
Jun/09/2016 18:19 Wrong answer on test 44 [main tests] → 18343827
Jun/09/2016 18:20 Time limit exceeded on test 50 [main tests] → 18343839
Jun/09/2016 18:21 Time limit exceeded on test 50 [main tests] → 18343863
Jun/09/2016 18:22 Time limit exceeded on test 50 [main tests] → 18343882
Jun/09/2016 18:24 Wrong answer on test 23 [main tests] → 18343912
Jun/09/2016 18:24 Wrong answer on test 40 [main tests] → 18343919
Jun/09/2016 18:33 Time limit exceeded on test 50 [main tests] → 18344067
Jun/09/2016 18:34 Time limit exceeded on test 50 [main tests] → 18344085
Jun/09/2016 18:34 Time limit exceeded on test 50 [main tests] → 18344097
Jun/09/2016 18:46 Time limit exceeded on test 50 [main tests] → 18344271
Jun/09/2016 18:50 Time limit exceeded on test 50 [main tests] → 18344326
Jun/09/2016 18:51 Time limit exceeded on test 50 [main tests] → 18344340
Jun/09/2016 19:06 Time limit exceeded on test 50 [main tests] → 18344571
Jun/10/2016 16:05 Wrong answer on test 10 [main tests] → 18358899
Jun/10/2016 16:06 Wrong answer on test 71 [main tests] → 18358925
Jun/10/2016 16:10 Wrong answer on test 71 [main tests] → 18358982
Jun/10/2016 16:12 Wrong answer on test 71 [main tests] → 18359017
Jun/10/2016 16:13 Wrong answer on test 71 [main tests] → 18359034
Jun/10/2016 16:14 Wrong answer on test 71 [main tests] → 18359049
Jun/10/2016 16:15 Time limit exceeded on test 1 [main tests] → 18359061
Jun/10/2016 16:17 Accepted [main tests] → 18359095

#include<stdio.h>
#include<vector>
#include<algorithm>
using namespace std;
#define SIZE 900
#define NS 167
int now[150400];
int flag[NS];
int dat[NS][150400];
int cnt[150400];
void resolve(int a)
{
	if (flag[a] != -1)
	{
		for (int i = 0; i < SIZE; i++)
		{
			dat[a][now[a*SIZE + i]] = 0;
			now[a*SIZE + i] = flag[a];
		}
		dat[a][flag[a]] = SIZE;
		flag[a] = -1;
	}
}
void update(int lb, int ub, int n)
{
	int b1 = lb / SIZE, b2 = ub / SIZE;
	resolve(b1);
	resolve(b2);
	if (b1 == b2)
	{
		for (int i = lb; i <= ub; i++)
		{
			dat[b1][now[i]]--;
			now[i] = n;
			dat[b1][n]++;
		}
	}
	else
	{
		for (int i = lb; i < (b1 + 1)*SIZE; i++)
		{
			dat[b1][now[i]]--;
			now[i] = n;
			dat[b1][n]++;
		}
		for (int i = b1 + 1; i < b2; i++)
		{
			flag[i] = n;
		}
		for (int i = b2*SIZE; i <= ub; i++)
		{
			dat[b2][now[i]]--;
			now[i] = n;
			dat[b2][n]++;
		}
	}
}
int getv(int p)
{
	if (flag[p / SIZE] != -1)return flag[p / SIZE];
	else return now[p];
}
int countn(int lb, int ub, int num)
{
	int b1 = lb / SIZE, b2 = ub / SIZE;
	int ret = 0;
	if (b1 == b2)
	{
		resolve(b1);
		for (int i = lb; i <= ub; i++)
		{
			ret += (now[i] == num);
		}
		return ret;
	}
	else
	{
		if (flag[b1] == -1)
		{
			if (lb%SIZE < SIZE / 2)
			{
				ret += dat[b1][num];
				for (int i = b1*SIZE; i < lb; i++)ret -= (now[i] == num);
			}
			else
			{
				for (int i = lb; i < (b1 + 1)*SIZE; i++)
				{
					ret += (now[i] == num);
				}
			}
		}
		else if (flag[b1] == num)ret += SIZE - (lb - b1*SIZE);
		for (int i = b1 + 1; i < b2; i++)
		{
			if (flag[i] != -1)ret += SIZE * (flag[i] == num);
			else ret += dat[i][num];
		}
		if (flag[b2] == -1)
		{
			if (ub%SIZE < SIZE / 2)
			{
				for (int i = b2*SIZE; i <= ub; i++)
				{
					ret += (now[i] == num);
				}
			}
			else
			{
				ret += dat[b2][num];
				for (int i = ub + 1; i < (b2 + 1)*SIZE; i++)ret -= (now[i] == num);
			}
		}
		else if (flag[b2] == num)ret += ub - b2*SIZE + 1;
		return ret;
	}
}
long long rrr = 2463534242LL;
int getrand(int mod)
{
	rrr = (rrr ^ (rrr << 13))&((1LL << 32) - 1); rrr = rrr ^ (rrr >> 17);
	rrr = (rrr ^ (rrr << 5))&((1LL << 32) - 1);
	return rrr%mod;
}
typedef pair<int, int>pii;
int main()
{
	int num, query, gen;
	scanf("%d%d%d", &num, &query, &gen);
	fill(now, now + 150400, -1);
	for (int i = 0; i < num; i++)
	{
		int z;
		scanf("%d", &z);
		z--;
		now[i] = z;
		dat[i / SIZE][z]++;
	}
	fill(flag, flag + NS, -1);
	for (int p = 0; p < query; p++)
	{
		int z, za, zb;
		scanf("%d%d%d", &z, &za, &zb);
		za--; zb--;
		if (z == 1)
		{
			int zc;
			scanf("%d", &zc);
			zc--;
			update(za, zb, zc);
		}
		else
		{
			int vec[100];
			int pt = 0;
			for (int i = 0; i < 100; i++)
			{
				int t = getv(getrand(zb - za + 1) + za);
				vec[pt++] = t;
				cnt[t]++;
			}
			vector<pii>v;
			for (int i = 0; i < 100; i++)
			{
				if (cnt[vec[i]] >= 3)v.push_back(make_pair(-cnt[vec[i]], vec[i]));
				cnt[vec[i]] = 0;
			}
			sort(v.begin(), v.end());
			vector<int>ans;
			for (int i = 0; i < min(int(v.size()),8); i++)
			{
				int t = countn(za, zb, v[i].second);
				if (t * 100 >= (zb - za + 1)*gen)
				{
					ans.push_back(v[i].second);
				}
			}
			printf("%d", ans.size());
			for (int i = 0; i < ans.size(); i++)printf(" %d", ans[i] + 1);
			printf("\n");
		}
	}
}

Codeforces 674F Bears and Juice

  • 問題概要

n人の人がいる。箱がたくさんあり、そのうち1個があたり。
毎日、各人は箱のうちいくつかを見て、その中に当たりが入っていた場合消える。p人まで消してよい。
i: 1~qについて、日数がi日のときに当たりの箱を特定できるような箱の個数の最大値を求めよ。

  • 解法

特定のk人が見た箱がx個あり、その中に当たりがあるとすると、k人が消えて当たりの箱の個数がx個になる。つまり、これは日数を1日、人数をk人減らし、箱の個数をx個にした状態と同じ。
よってDP[消えた人数][残り日数]: それ以降で特定できる箱の個数の最大値x としてDPができる。

このDPはO(p^2q)で間に合わないが、DPの遷移は残り日数によらず、次の日のDP配列の線形結合なので、同じ行列の累乗の形でかける。

さらにこの行列は上三角(消えた人の人数は減らないため)なので、求める答えはあるp次多項式に1~qを代入した値となる。
この多項式の小さいところでの値は普通にDPをすると求まるので、補間をする。この場合隣同士の差をp回とってやると定数関数になるので、imos法っぽくやるとO(pq)で全部の値が求められる。

全体でO(p^3+pq)となる。ゆっくり書いたらバグらなかったのでOK。DPしようと思うのが難しいと思う。(解説を昔5秒眺めたことがあってDPという情報だけ知っていた)

#include<stdio.h>
#include<vector>
#include<algorithm>
using namespace std;
typedef long long ll;
typedef pair<ll, ll>pii;
pii num2(ll a)
{
	int ret = 0;
	for (;;)
	{
		if (a % 2 != 0)return make_pair(ret, a);
		a /= 2;
		ret++;
	}
}
ll mat[131][131];
ll rev[200];
ll now[131];
ll imos[131][131];
int main()
{
	ll num, gen, query;
	scanf("%I64d%I64d%I64d", &num, &gen, &query);
	gen = min(gen, num - 1);
	for (int i = 1; i < 200; i++)
	{
		if (i % 2 != 0)
		{
			for (ll j = 0;; j++)
			{
				if (((1LL << 32)*j + 1) % i == 0)
				{
					rev[i] = ((1LL << 32)*j + 1) / i;
					break;
				}
			}
		}
	}
	for (int i = 0; i <= gen; i++)
	{
		for (int j = i; j <= gen; j++)
		{
			ll n = num - i;
			ll p = j - i;
			int c = 0;
			ll now = 1;
			for (int k = 1; k <= p; k++)
			{
				pii z1 = num2(n - k + 1);
				pii z2 = num2(k);
				c += z1.first - z2.first;
				now *= z1.second*rev[z2.second] & ((1LL << 32) - 1);
				now &= ((1LL << 32) - 1);
			}
			mat[i][j] = now << (c % 32);
		}
	}
	fill(now, now + 131, 1);
	for (int i = 0; i <= gen; i++)
	{
		int n[131];
		fill(n, n + 131, 0);
		for (int j = 0; j <= gen; j++)
		{
			for (int k = 0; k <= gen; k++)
			{
				n[j] += mat[j][k] * now[k];
				n[j] &= ((1LL << 32) - 1);
			}
		}
		imos[0][i] = now[0];
		for (int j = 0; j <= gen; j++)now[j] = n[j];
	}
	for (int i = 0; i < gen; i++)
	{
		for (int j = 0; j < gen - i; j++)
		{
			imos[i + 1][j] = (imos[i][j + 1] - imos[i][j] + (1LL << 32))&((1LL << 32) - 1);
		}
	}
	ll t = imos[gen][0];
	ll ans[131];
	for (int i = 0; i <= gen; i++)ans[i] = imos[i][0];
	ll ret = 0;
	for (int i = 1; i <= query; i++)
	{
		ll n[131];
		fill(n, n + 131, 0);
		n[gen] = t;
		for (int j = gen - 1; j >= 0; j--)
		{
			n[j] = (ans[j + 1] + ans[j])&((1LL << 32) - 1);
		}
		for (int j = 0; j <= gen; j++)ans[j] = n[j];
		ret ^= (ans[0] * i)&((1LL << 32) - 1);
	}
	printf("%I64d\n", ret);
}

Codeforces 666E Forensic Examination

  • 問題概要:

文字列Sと、文字列X_1,...,X_Nが与えられる。次のクエリをQ個処理せよ。

クエリ: l,r,p,qが与えられるので、X_lからX_rのなかで文字列[s_p...s_q]が部分文字列として一番多く現れるものとその回数を求めよ。

  • 解法:

suffix array+LCP+sparse tableでクエリを区間の最頻値を求める形にしておく。この区間たちはラミナー(交わるなら一方が他方に含まれる)なので、マージテクをしながらクエリを区間が小さい順に見ていくとsegtreeの操作で解ける。

サイズがでかいので前半パートにlogを2つつけると死にます。後半パートは50000につくので大丈夫。

こういう実装遅すぎるし演習量が足りません。精進します

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<string>
#include<iostream>
#include<set>
#include<time.h>
using namespace std;
vector<int>sais(vector<int>vec)//1-originにして末尾に0を加えてから呼ぶのを忘れないこと
{
	if (vec.empty())
	{
		vector<int>v;
		return v;
	}
	vector<int>dat;
	dat.resize(vec.size());
	dat[vec.size() - 1] = 0;
	for (int i = vec.size() - 2; i >= 0; i--)
	{
		if (vec[i]>vec[i + 1])dat[i] = 1;
		else if (vec[i]<vec[i + 1])dat[i] = 0;
		else dat[i] = dat[i + 1];
	}
	vector<vector<int> >sa;
	int maxi = 0;
	for (int i = 0; i<vec.size(); i++)maxi = max(maxi, vec[i]);
	sa.resize(maxi + 1);
	vector<int>vnum;
	vnum.resize(maxi + 1);
	for (int i = 0; i<vec.size(); i++)vnum[vec[i]]++;
	vector<bool>islms;
	islms.resize(vec.size());
	fill(islms.begin(), islms.end(), false);
	for (int i = 0; i <= maxi; i++)
	{
		sa[i].resize(vnum[i]);
		fill(sa[i].begin(), sa[i].end(), -1);
	}
	vector<int>pt1, pt2;
	pt1.resize(maxi + 1);
	pt2.resize(maxi + 1);
	for (int i = 0; i <= maxi; i++)
	{
		pt2[i] = vnum[i] - 1;
	}
	for (int i = vec.size() - 1; i >= 1; i--)
	{
		if ((dat[i - 1] == 1 && dat[i] == 0) || i == vec.size() - 1)
		{
			sa[vec[i]][pt2[vec[i]]--] = i;
			islms[i] = true;
		}
	}
	for (int i = 0; i <= maxi; i++)
	{
		for (int j = 0; j<vnum[i]; j++)
		{
			if (sa[i][j]>0)
			{
				if (dat[sa[i][j] - 1] == 1)
				{
					sa[vec[sa[i][j] - 1]][pt1[vec[sa[i][j] - 1]]++] = sa[i][j] - 1;
				}
			}
		}
	}
	for (int i = 1; i <= maxi; i++)
	{
		for (int j = pt2[i] + 1; j<vnum[i]; j++)sa[i][j] = -1;
		pt2[i] = vnum[i] - 1;
	}
	for (int i = maxi; i >= 0; i--)
	{
		for (int j = vnum[i] - 1; j >= 0; j--)
		{
			if (sa[i][j]>0)
			{
				if (dat[sa[i][j] - 1] == 0)
				{
					sa[vec[sa[i][j] - 1]][pt2[vec[sa[i][j] - 1]]--] = sa[i][j] - 1;
				}
			}
		}
	}
	vector<int>d;
	d.resize(vec.size());
	int cnt = 0;
	vector<int>bef;
	bef.push_back(0);
	bool fl = false;
	for (int i = 0; i <= maxi; i++)
	{
		for (int j = 0; j<vnum[i]; j++)
		{
			if (islms[sa[i][j]])
			{
				vector<int>zk;
				zk.push_back(vec[sa[i][j]]);
				for (int k = sa[i][j] + 1; k<vec.size(); k++)
				{
					zk.push_back(vec[k]);
					if (islms[k])break;
				}
				if (zk != bef)cnt++;
				else if (vec[sa[i][j]] != 0)fl = true;
				d[sa[i][j]] = cnt;
				bef = zk;
			}
		}
	}
	vector<int>vt;
	for (int i = 0; i<vec.size(); i++)if (islms[i])vt.push_back(d[i]);
	vector<int>gv;
	vector<int>nv;
	if (fl)
	{
		gv = sais(vt);
		vector<int>v;
		for (int i = 0; i<vec.size(); i++)
		{
			if (islms[i])v.push_back(i);
		}
		for (int i = 0; i<gv.size(); i++)
		{
			nv.push_back(v[gv[i]]);
		}
	}
	else
	{
		gv = vt;
		nv.resize(gv.size());
		int pt = 0;
		for (int i = 0; i<gv.size(); i++)
		{
			for (;;)
			{
				if (islms[pt])
				{
					nv[gv[i]] = pt;
					pt++;
					break;
				}
				pt++;
			}
		}
	}
	for (int i = 0; i <= maxi; i++)
	{
		fill(sa[i].begin(), sa[i].end(), -1);
		pt1[i] = 0;
		pt2[i] = vnum[i] - 1;
	}
	for (int i = nv.size() - 1; i >= 0; i--)
	{
		sa[vec[nv[i]]][pt2[vec[nv[i]]]--] = nv[i];
	}
	for (int i = 0; i <= maxi; i++)
	{
		for (int j = 0; j<vnum[i]; j++)
		{
			if (sa[i][j]>0)
			{
				if (dat[sa[i][j] - 1] == 1)
				{
					sa[vec[sa[i][j] - 1]][pt1[vec[sa[i][j] - 1]]++] = sa[i][j] - 1;
				}
			}
		}
	}
	for (int i = 1; i <= maxi; i++)
	{
		for (int j = pt2[i] + 1; j<vnum[i]; j++)sa[i][j] = -1;
		pt2[i] = vnum[i] - 1;
	}
	for (int i = maxi; i >= 0; i--)
	{
		for (int j = vnum[i] - 1; j >= 0; j--)
		{
			if (sa[i][j]>0)
			{
				if (dat[sa[i][j] - 1] == 0)
				{
					sa[vec[sa[i][j] - 1]][pt2[vec[sa[i][j] - 1]]--] = sa[i][j] - 1;
				}
			}
		}
	}
	vector<int>ret;
	for (int i = 0; i <= maxi; i++)
	{
		for (int j = 0; j<vnum[i]; j++)
		{
			ret.push_back(sa[i][j]);
		}
	}
	return ret;
}
vector<int>calclcp(vector<int>str, vector<int>sa)//lcp: SAのi-1番目とi番目のlcp
{
	vector<int>rsa;
	rsa.resize(sa.size());
	for (int i = 0; i < sa.size(); i++)rsa[sa[i]] = i;
	vector<int>lcp;
	lcp.resize(sa.size());
	int now = 1;
	for (int i = 0; i < str.size() - 1; i++)
	{
		if (now != 0)now--;
		for (;;)
		{
			if (str[i + now] == str[sa[rsa[i] - 1] + now])now++;
			else
			{
				lcp[rsa[i]] = now;
				break;
			}
		}
	}
	return lcp;
}
#define SIZE 1048576
class segtree
{
public:
	int seg[SIZE * 2];
	void init()
	{
		for (int i = 1; i < SIZE * 2; i++)seg[i] = 1000000000;
	}
	void update(int a, int b)
	{
		a += SIZE;
		seg[a] = min(seg[a], b);
		for (;;)
		{
			a /= 2;
			if (a == 0)break;
			seg[a] = min(seg[a * 2], seg[a * 2 + 1]);
		}
	}
	int get(int beg, int end, int node, int lb, int ub)
	{
		if (ub < beg || end < lb)return 1000000000;
		if (beg <= lb&&ub <= end)return seg[node];
		return min(get(beg, end, node * 2, lb, (lb + ub) / 2), get(beg, end, node * 2 + 1, (lb + ub) / 2 + 1, ub));
	}
};
class sparsetable
{
public:
	int seg[20][SIZE];
	int dat[SIZE];
	void init(vector<int>vec)
	{
		for (int i = 1; i < SIZE; i++)
		{
			int c = 0, now = i;
			for (;;)
			{
				if (now == 0)break;
				now /= 2;
				c++;
			}
			dat[i] = c - 1;
		}
		for (int i = 0; i < vec.size(); i++)seg[0][i] = vec[i];
		for (int i = vec.size(); i < SIZE; i++)seg[0][i] = 1000000000;
		for (int i = 1; i < 20; i++)
		{
			for (int j = 0; j < SIZE; j++)
			{
				seg[i][j] = min(seg[i - 1][j], seg[i - 1][j + (1 << (i - 1))]);
			}
		}
	}
	int get(int lb, int ub)
	{
		int t = ub - lb + 1;
		if (t <= 0)return 0;
		return min(seg[dat[t]][lb], seg[dat[t]][ub - (1 << dat[t]) + 1]);
	}
};
sparsetable sp;
//segtree tree;
vector<int>str, sa, lcp, rsa;
void init(string s)
{
	for (int i = 0; i < s.size(); i++)str.push_back(s[i] - 'a' + 1);
	str.push_back(0);
	sa = sais(str);
	lcp = calclcp(str, sa);
	//tree.init();
	sp.init(lcp);
	rsa.resize(sa.size());
	for (int i = 0; i < sa.size(); i++)rsa[sa[i]] = i;
	//for (int i = 0; i < str.size(); i++)tree.update(i, lcp[i]);
}
int getlcp(int a, int b)//a文字目からとb文字目からのlcpの長さ
{
	return sp.get(min(rsa[a], rsa[b]) + 1, max(rsa[a], rsa[b]));
}
typedef pair<int, int>pii;
typedef pair<pii, int>pi3;
typedef pair<pi3, pii>pi5;
pii ans[500000];
class segtreedy
{
public:
	vector<int>seg;
	vector<int>idx;
	vector<int>lko;
	vector<int>rko;
	vector<pii>vec;
	void init()
	{
		seg.push_back(0);
		lko.push_back(-1);
		rko.push_back(-1);
		idx.push_back(-1);
	}
	void add(int pl, int node, int lb, int ub, int num)
	{
		if (ub < pl || pl < lb)return;
		if (lb == ub)
		{
			seg[node] += num;
			idx[node] = -pl;
			vec.push_back(make_pair(pl, num));
			return;
		}
		if (lko[node] == -1)
		{
			lko[node] = seg.size();
			seg.push_back(0);
			lko.push_back(-1);
			rko.push_back(-1);
			idx.push_back(-1);
		}
		if (rko[node] == -1)
		{
			rko[node] = seg.size();
			seg.push_back(0);
			lko.push_back(-1);
			rko.push_back(-1);
			idx.push_back(-1);
		}
		add(pl, lko[node], lb, (lb + ub) / 2, num);
		add(pl, rko[node], (lb + ub) / 2 + 1, ub, num);
		seg[node] = max(seg[lko[node]], seg[rko[node]]);
		if (seg[node] == seg[lko[node]])idx[node] = idx[lko[node]];
		else idx[node] = idx[rko[node]];
	}
	pii get(int beg, int end, int node, int lb, int ub)
	{
		if (ub < beg || end < lb)return make_pair(0, -1);
		if (beg <= lb&&ub <= end)
		{
			//if (idx[node] == -1)return make_pair(seg[node], -lb);
			return make_pair(seg[node], idx[node]);
		}
		if (lko[node] == -1)
		{
			lko[node] = seg.size();
			seg.push_back(0);
			lko.push_back(-1);
			rko.push_back(-1);
			idx.push_back(-1);
		}
		if (rko[node] == -1)
		{
			rko[node] = seg.size();
			seg.push_back(0);
			lko.push_back(-1);
			rko.push_back(-1);
			idx.push_back(-1);
		}
		return max(get(beg, end, lko[node], lb, (lb + ub) / 2), get(beg, end, rko[node], (lb + ub) / 2 + 1, ub));
	}
	void destroy()
	{
		seg.clear();
		lko.clear();
		rko.clear();
		vec.clear();
	}
};
segtreedy dat[600000];
int main()
{
	time_t clll = clock();
	string s;
	cin >> s;
	vector<string>vec;
	string conc;
	int num;
	scanf("%d", &num);
	vector<int>toidx;
	for (int i = 0; i < num; i++)
	{
		string z;
		cin >> z;
		vec.push_back(z);
		conc += z;
		conc.push_back('z' + 1);
		for (int j = 0; j < z.size(); j++)toidx.push_back(i + 1);
		toidx.push_back(0);
	}
	int spt = conc.size();
	conc += s;
	for (int i = 0; i < s.size(); i++)toidx.push_back(0);
	init(conc);
	int query;
	scanf("%d", &query);
	vector<pi5>que;
	for (int i = 0; i < query; i++)
	{
		int za, zb, zc, zd;
		scanf("%d%d%d%d", &za, &zb, &zc, &zd);
		zc--, zd--;
		int lb, ub;
		int beg = 0, end = rsa[spt + zc];
		for (;;)
		{
			if (beg == end)break;
			int med = (beg + end) / 2;
			if (getlcp(sa[med], spt + zc) >= zd - zc + 1)end = med;
			else beg = med+1;
		}
		lb = beg;
		beg = rsa[spt + zc], end = conc.size() - 1;
		for (;;)
		{
			if (beg == end)break;
			int med = (beg + end+1) / 2;
			if (getlcp(sa[med], spt + zc) >= zd - zc + 1)beg = med;
			else end = med-1;
		}
		ub = beg;
		que.push_back(make_pair(make_pair(make_pair(ub - lb + 1, lb), i), make_pair(za, zb)));
	}
	//for (int i = 0; i < conc.size(); i++)printf("%d", i % 10); printf("\n");
	//cout << conc << endl;
	//for (int i = 0; i < conc.size(); i++)printf("%d %c %d %d\n", i, conc[sa[i]], sa[i], lcp[i]);
	sort(que.begin(), que.end());
	for (int i = 0; i < conc.size(); i++)
	{
		dat[i].init();
		if (toidx[sa[i]] != 0)dat[i].add(toidx[sa[i]], 0, 0, SIZE - 1, 1);
	}
	//if (double(clock() - clll) / CLOCKS_PER_SEC>2.9)return 0;
	//for (int i = 0; i < conc.size(); i++)printf("%d ", toidx[sa[i]]);
	//printf("\n");
	set<pii>se;
	for (int i = 0; i < conc.size(); i++)if (toidx[sa[i]] != 0)se.insert(make_pair(i, i));
	se.insert(make_pair(conc.size(), conc.size()));
	se.insert(make_pair(conc.size() + 1, conc.size()));
	for (int i = 0; i < que.size(); i++)
	{
		int now = que[i].first.first.second;
		int end = que[i].first.first.second + que[i].first.first.first - 1;
		//printf("%d %d\n", now, end);
		set<pii>::iterator it = se.lower_bound(make_pair(now,-1));
		for (;;)
		{
			int pt1 = (*it).second;
			it++;
			if ((*it).first >= end + 1)break;
			int pt2 = (*it).second;
			se.erase(it++);
			it--;
			se.erase(it++);
			if (dat[pt1].vec.size() > dat[pt2].vec.size())swap(pt1, pt2);
			for (int j = 0; j < dat[pt1].vec.size(); j++)
			{
				//printf(" %d %d\n", dat[pt1].vec[j].first, dat[pt1].vec[j].second);
				dat[pt2].add(dat[pt1].vec[j].first, 0, 0, SIZE - 1, dat[pt1].vec[j].second);
			}
			dat[pt1].destroy();
			se.insert(make_pair(now, pt2));
			it--;
		}
		it = se.lower_bound(make_pair(now, -1));
		if ((*it).first > end)
		{
			ans[que[i].first.second] = make_pair(0, -que[i].second.first);
			continue;
		}
		pii t = dat[(*it).second].get(que[i].second.first, que[i].second.second, 0, 0, SIZE - 1);
		if (t.first == 0)t.second = -que[i].second.first;
		ans[que[i].first.second] = t;
		//for (int i = 0; i < 5; i++)printf("%d %d  ", dat[(*it).second].get(i, i, 0, 0, SIZE - 1).first, dat[(*it).second].get(i, i, 0, 0, SIZE - 1).second); printf("\n");
	}
	for (int i = 0; i < query; i++)
	{
		printf("%d %d\n", -ans[i].second, ans[i].first);
	}
}

CF 296E Triangles 3000

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

問題概要: 平面上に直線が3000本くらいあるので、その中の3直線をランダムに選んだ時の囲まれる面積の期待値を求めよ。


解法: 三角形ABCの面積は、原点をOとすると三角形OAB,OBC,OCAの面積の符号付和で書ける。OABにつく符号は、Cが直線ABに対してOと同じ側にあるかどうかで決まる。OBC,OCAも同様。

直線を1本固定して、その直線を使う三角形について、その直線上の2点とOからなる三角形の面積の符号付和の総和を求めるというのをすべての直線についてやればいい。
これは、直線を固定すれば、その直線上の2点とOからなる三角形の高さは常に等しいので底辺の総和を求めればよく、その直線とほかの直線との交点とその交わりの偏角を全列挙して偏角順にソートし、偏角の小さい順に点を追加しながら、適切なデータ構造で座標の和を扱えばいい。全体でO(N^2log N)となり、通る。実装では、直線を扱うごとにそれが水平になるように図全体を回転すると楽。

#include<stdio.h>
#include<vector>
#include<algorithm>
#include<math.h>
#include<stdlib.h>
using namespace std;
typedef pair<double, double>pdd;
typedef pair<pdd, double>pd3;
#define SIZE 4096
class BIT
{
public:
	double bit[SIZE + 1];
	void init()
	{
		fill(bit, bit + SIZE + 1, 0.0);
	}
	void add(int a, double b)
	{
		a++;
		for (;;)
		{
			if (a > SIZE)return;
			bit[a] += b;
			a += a&-a;
		}
	}
	double get(int a)
	{
		a++;
		double ret = 0.0;
		for (;;)
		{
			if (a == 0)return ret;
			ret += bit[a];
			a -= a&-a;
		}
	}
};
BIT bi1, bi2;
double pi = 3.1415926535897932384626433832;
double calc(vector<pd3>dat, pd3 lin)
{
	double a = lin.first.first, b = lin.first.second, c = lin.second;
	double vx = a / sqrt(a*a + b*b), vy = b / sqrt(a*a + b*b);
	vector<pdd>vec;
	vector<double>zat;
	//printf("%lf %lf\n", vx, vy);
	for (int i = 0; i < dat.size(); i++)
	{
		double za = dat[i].first.first, zb = dat[i].first.second, zc = dat[i].second;
		double arg = atan2(a, b) - atan2(za, zb);
		if (arg > pi)arg -= pi;
		if (arg > pi)arg -= pi;
		if (arg < 0)arg += pi;
		if (arg < 0)arg += pi;
		double x = (zb*c - b*zc) / (a*zb - b*za), y = (c*za - a*zc) / (b*za - a*zb);
		//printf(" %lf %lf %lf\n", x, y, arg);
		vec.push_back(make_pair(arg, x*vy - y*vx));
		zat.push_back(x*vy - y*vx);
	}
	sort(vec.begin(), vec.end());
	sort(zat.begin(), zat.end());
	bi1.init();
	bi2.init();
	double d = c / sqrt(a*a + b*b);
	double ret = 0;
	for (int i = 0; i < dat.size(); i++)
	{
		//printf("  %lf %lf\n", vec[i].first, vec[i].second);
		int low = lower_bound(zat.begin(), zat.end(), vec[i].second) - zat.begin();
		ret -= bi2.get(low)*vec[i].second - bi1.get(low);
		ret += (bi1.get(SIZE - 1) - bi1.get(low)) - (bi2.get(SIZE - 1) - bi2.get(low))*vec[i].second;
		bi1.add(low, vec[i].second);
		bi2.add(low, 1.0);
	}
	//printf("%lf %lf\n\n", ret, d);
	return ret*d / 2.0;
}
int main()
{
	int num;
	scanf("%d", &num);
	vector<pd3>vec;
	for (int i = 0; i < num; i++)
	{
		double za, zb, zc;
		scanf("%lf%lf%lf", &za, &zb, &zc);
		vec.push_back(make_pair(make_pair(za, zb), zc));
	}
	double ans = 0.0;
	for (int i = 0; i < num; i++)
	{
		vector<pd3>z;
		for (int j = 0; j < num; j++)
		{
			if (j != i)z.push_back(vec[j]);
		}
		ans += calc(z, vec[i]);
	}
	ans *= 6.0;
	ans /= num;
	ans /= (num - 1);
	ans /= (num - 2);
	printf("%lf\n", ans);
}