POI: Ploughing

https://szkopul.edu.pl/problemset/problem/6YiP6JA5U15hY94pLwuHoYPg/site/?key=statement

実践的アルゴリズム貪欲編です。

少し観察すると(列をstripする回数==M)または(行をstripする回数==N)が成立することがわかります。
対称なので(列をstripする回数==M)とします。
行について上からstripする回数と下からstripする回数を固定すれば、列は貪欲的にstripを行えます。
なのでdp[上からstripする回数][下からstripする回数]=make_pair(左から何回stripできるか,右から何回stripできるか) を埋めればよくて、それはならし計算量O((N+M)^2)で行うことができます。

ほそながいところより難しいということになっていましたが、こっちの方がまだ考えやすいと思いました。

int K, N, M;
int XS[2010][2010];//x-axis way sum,[ )
int YS[2010][2010];//y-axis way sum,[ )
pi dp[2010][2010];//[0, l) and [r, N)

pi narrow(int xl, int xr, int yl, int yr) { //xl, xr:constant, narrow yl and yr. [0, yl) and [yr, M)
	while(yr - yl >= 1 && XS[xr][yl] - XS[xl][yl] <= K) {
		yl++;
	}
	while(yr - yl >= 1 && XS[xr][yr - 1] - XS[xl][yr - 1] <= K) {
		yr--;
	}
	return pi(yl, yr);
}

int solve_sub() {
	int res = inf;
	rep(i, 0, N + 1) 
		rep(j, i + 1, N + 1) dp[i][j] = pinf;
	dp[0][N] = narrow(0, N, 0, M);
	rep(i, 0, N + 1) {
		rer(k, N + 1, 2) {
			int j = i + k;
			if(j > N || dp[i][j] == pinf) continue;
			int yl = dp[i][j].fst, yr = dp[i][j].sec;
			if(dp[i + 1][j] == pinf && YS[i][yr] - YS[i][yl] <= K) {
				dp[i + 1][j] = narrow(i + 1, j, yl, yr);
			}
			if(dp[i][j - 1] == pinf && YS[j - 1][yr] - YS[j - 1][yl] <= K) {
				dp[i][j - 1] = narrow(i, j - 1, yl, yr);
			}
		}
	}
	rep(i, 0, N + 1) {
		rep(j, 0, N + 1) {
			if(dp[i][j] == pinf) continue;
			if(dp[i][j].sec - dp[i][j].fst == 0) {
				MIN(res, M + i + N - j);
			}
		}
	}
	return res;
}

void solve() {
	cin >> K >> M >> N;
	rep(i, 0, N) {
		rep(j, 0, M) cin >> XS[i + 1][j];
	}
	rep(i, 0, N) 
		rep(j, 0, M) YS[i][j + 1] = YS[i][j] + XS[i + 1][j];
	rep(j, 0, M)
		rep(i, 1, N + 1) XS[i + 1][j] += XS[i][j];
	int res = inf;
	MIN(res, solve_sub());
	rep(i, 0, max(N, M) + 1) {
		rep(j, i + 1, max(N, M) + 1) {
			swap(XS[i][j], XS[j][i]);
			swap(YS[i][j], YS[j][i]);
		}
	}
	swap(XS, YS);
	swap(N, M);
	MIN(res, solve_sub());
	cout << res << "\n";
}

無駄にMLEが厳しいのやめて…

POJ 1284: Primitive Roots

http://poj.org/problem?id=1284

答えはφ(p-1)個になります。証明しましょう。

準備.n(|p-1)についてx^n≡1(modp)はちょうどn個の解を持つ。

証明.p-1=nkとしたとき、
x^(p-1)-1=(x^n-1)(x^(n(k-1))+x^(n(k-1)-1)+...+1)であり、
x^(p-1)-1≡0(modp)の解はp-1=nk個。
x^n-1≡0(modp)の解はn個以下
x^(n(k-1))+x^(n(k-1)-1)+...+1≡0の解はn(k-1)個以下なので
x^n-1≡(modp)の解はn個でなければならない。

n個の解x1,...,xnについてx1^a≡1(modp)を満たす最小の自然数aを位数と定義すると、位数はnの約数しかありえないので、
位数aの解の個数をF(a)と置くと、Σ(a|n)F(a)=nとなる。これがp-1の任意の約数nについて成立している。
またF(1)=1
F(p)=p-1
F(p^2)=p^2-p(pは素数)
F(ab)=F(a)F(b)(aとbは互いに素)が容易に示せ、帰納的にF(a)=φ(a)が成立する。
よって位数p-1の解の個数はF(p-1)=φ(p-1)となる。


x^k≡1(modp)となる集合A(k)を包除原理して求めてもよい。φ(n)を求めるときと合同な集合が出てくる。

ちなみに任意の自然数についてはφ(φ(N))個となります。


ほかにもこの証明をする過程でいろいろ考えた定理を挙げておきます。

  • AとNが互いに素の時、A^φ(N)≡1(modN)(証明はA*2A*...*φ(N)A≡1*2*...*φ(N)(modN)より)

(AとNが互いに素ではないときはgcd(A,N)に含まれている素因数を除いたA',N'について適応する)

  • A|Bの時φ(A)|φ(B) (φ(A)=A*(p1-1)/p1*(p2-1)/p2から簡単に示せる)
  • 十分大きなKをとればA^(K+φ(N))≡A^K(modN)(フェルマーの定理と中国剰余定理から)

これよりA^K(modN)の値はK modφ(N)に依存。

  • 素数Pの原始根の1つをgとしたとき1<=z<=p-1なる自然数zについてz≡g^i(modP)となるiが存在する。(原始根の定義そのもの)

これを利用して原始根の数がφ(φ(N))となることを証明する方法もある。

  • 原始根を効率よく見つけるアルゴリズムはない。(だから暗号とかがうまくいく)
int N;
int phi[100010];

int main() {
	rep(i, 1, 100000) phi[i] = i;
	rep(i, 2, 100000) {
		if(phi[i] != i) continue;
		for(int j = i; j < 100000; j += i) {
			phi[j] -= phi[j] / i;
		}
	}
	while(scanf("%d", &N) != EOF) {
		printf("%d\n", phi[N - 1]);
	}
}

POJ 2115: C Looooopsなど

POJ 2115:
http://poj.org/problem?id=2115

A+Cx=B+(2^K)yの解(x,y)のうちxが非負で最小のものを求めればよい。これはユークリッドの互除法やるだけで求められるが、少し遠回りに書いたらREを起こしてえぇ…ってなりました。たぶんオーバーフロー起こしてましたね。シンプルに書いてAC。

ll extgcd(ll a, ll b, ll& x, ll& y) {
	if(b == 0) {
		x = 1; y = 0;
		return a;
	}
	ll res = extgcd(b, a % b, x, y);
	ll tx = x, ty = y;
	x = -ty; y = -tx - (a / b) * ty;
	return res;
}

int main() {
	while(true) {
		ll A, B, C, K;
		scanf("%lld%lld%lld%lld", &A, &B, &C, &K);
		if(A == 0 && B == 0 && C == 0 && K == 0) break;
		K = (1LL << K);
		ll x, y;
		ll G = extgcd(K, C, y, x);
		if((A - B) % G != 0) {
			printf("FOREVER\n");
			continue;
		}
		K /= G; C /= G;
		x = x * ((A - B) / G);
		printf("%lld\n", (x % K + K) % K);
	}
	return 0;
}

POJ 2345:
http://poj.org/problem?id=2345

いっけんムリそうですが、行列の形で書き下すと連立方程式解くだけでよいことがわかるので掃き出し法などすればよいです。
ムリそうだなぁと思ったらフローとか数学とかで考えるのがよさそう。

bool gauss_jordan(mat& A, vi& B) {
	int n = sz(B);
	rep(j, 0, n) {
		int at = -1;
		rep(i, j, n) {
			if(A[i][j] == 1) {
				at = i;
				break;
			}
		}
		if(at == -1) continue;
		swap(A[j], A[at]); swap(B[j], B[at]);
		rep(i, 0, n) {
			if(i == j || A[i][j] == 0) continue;
			B[i] = B[i] ^ B[j];
			rep(k, 0, n) {
				A[i][k] = A[i][k] ^ A[j][k];
			}
		}
	}
	rep(i, 0, n) {
		bool found = false;
		rep(j, 0, n) {
			if(A[i][j]) {
				found = true; break;
			}
		}
		if(!found && B[i]) return false;
	}
	return true;
}

int N;

int main() {
	scanf("%d", &N);
	mat A(N, vi(N, 0));
	rep(i, 0, N) {
		while(true) {
			int a;
			scanf("%d", &a); a--;
			if(a == -2) break;
			A[a][i] = 1;
		}
	}
	vi B(N, 1);
	if(gauss_jordan(A, B)) {
		vi ans;
		rep(i, 0, N) {
			if(B[i]) ans.pb(i + 1);
		}
		rep(i, 0, sz(ans)) {
			printf("%d%c", ans[i], i != sz(ans) - 1 ? ' ' : '\n');
		}
	}
	else printf("No solution\n");
	return 0;
}

POJ 1150: The Last Non-zero Digit

http://poj.org/submit?problem_id=1150

初中国剰余定理かもしれない。
(X/10^k)mod10を求めればいいです。ここでkはXが10で割り切れる回数です。
中国剰余定理より
(X/10^k)mod10 <=> (X/10^k)mod2, (X/10^k)mod5となります。
対称性よりmod5についてのみ注目して考えましょう。
X=a*5^eとします。ここでa mod5 != 0です。a mod5とeは蟻本の方法で求められます。するとe>=kが成立するのは明らかです。
e>kのとき、(X/10^k)mod5 = 0です。
e=kのとき、(X/10^k)≡(X/5^k)*(1/2^k)≡a*inv[2^k](mod5)となるので求められます。

int fact[6][11];
int inv[6][11];

ll mod_pow(ll a, ll n, ll mod) {
	if(n == 0) return 1;
	ll res = mod_pow(a, n / 2, mod);
	if(n % 2 == 0) return res * res % mod;
	else return a * res % mod * res % mod;
}


int mod_fact(int n, int p, int &e) {
	e = 0;
	if(n == 0) return 1;
	int res = mod_fact(n / p, p, e);
	e += n / p;

	if((n / p) % 2 != 0) return res * (p - fact[p][n % p]) % p;
	else return res * fact[p][n % p] % p;
}

int N, M;
int ans[5][10];

int main() {
	fact[2][0] = fact[2][1] = 1;
	inv[2][1] = 1;
	fact[5][0] = fact[5][1] = 1;
	inv[5][1] = 1;
	rep(i, 2, 5) {
		fact[5][i] = fact[5][i - 1] * i % 5;
		inv[5][i] = 5 - inv[5][5 % i] * (5 / i) % 5;
	}
	rep(i, 0, 10) {
		ans[i % 2][i % 5] = i;
	}
	while(scanf("%d%d", &N, &M) != EOF) {
		int e2_1, e2_2;
		int e5_1, e5_2;
		int a2_1 = mod_fact(N, 2, e2_1);
		int a2_2 = mod_fact(N - M, 2, e2_2);
		int a5_1 = mod_fact(N, 5, e5_1);
		int a5_2 = mod_fact(N - M, 5, e5_2);
		int a2 = a2_1 * inv[2][a2_2] % 5;
		int a5 = a5_1 * inv[5][a5_2] % 5;
		int e2 = e2_1 - e2_2;
		int e5 = e5_1 - e5_2;
		if(e2 > e5) {
			a5 = a5 * inv[5][mod_pow(2, e5, 5) % 5] % 5;
			printf("%d\n", ans[0][a5]);
		}
		else if(e2 == e5) {
			a5 = a5 * inv[5][mod_pow(2, e5, 5) % 5] % 5;
			a2 = a2 * inv[2][mod_pow(5, e2, 2) % 2] % 2;
			printf("%d\n", ans[a2][a5]);
		}
		else {
			a2 = a2 * inv[2][mod_pow(5, e2, 2) % 2] % 2;
			printf("%d\n", ans[a2][0]);
		}
	}
}

CODE FESTIVAL 2017 qual B,E: Popping Balls

http://code-festival-2017-qualb.contest.atcoder.jp/tasks/code_festival_2017_qualb_e

赤のボールがA,青のボールがB個の状態から選ぶ方法の数を再帰的に求めます。これを(A, B)と定義しておきましょう。
赤を選ぶ場合単に(A-1,B)とすればいいです。青を選ぶ場合t=B+1にします。そうするとボールの数がt個になるまで赤と青のボールを好きなように選べます。それで赤がp個、青がq個(p+q=t+1)になったとしましょう。赤を使った場合は(p-1,q)、青を使った場合はs=q+1としてまたボールがq+1個になるまで好きなように選べます。q+1個になったらあとは赤を全部取ってから青をとるしかないです。なのでp,qをすべて試して(p,q)と(t+1,Bからp,qへの経路数)の積の和を求めればいいです。後者は単なるコンビネーションで、前者はdpで前処理計算できます。

区間の変な問題に言い換えてしまい惨敗でした。コンビネーションがたくさん出そうだということはわかったのですが、経路数に言い換えることはできませんでした。方針の変更は難しい。

あと再帰的に考えることの重要性も学びました。自分はdpを配る形でよく書くので再帰的に処理していることを忘れがちです。

いい言い換え+再帰的思考ができれば解けた問題でした。

int A, B;
ll C[2020][2020];
ll CS[2020][2020];
ll dp[2020][2020];

void C_init(int n) {
	n++;
	rep(i, 0, n) C[i][0] = 1;
	rep(i, 1, n) {
		rep(j, 1, i + 1) {
			C[i][j] = C[i - 1][j] + C[i - 1][j - 1];
			C[i][j] %= mod;
		}
	}
}

void solve() {
	cin >> A >> B;
	C_init(2010);
	rep(i, 0, 2010) {
		rep(j, 0, 2010) {
			CS[i][j + 1] = CS[i][j] + C[i][j];
			CS[i][j + 1] %= mod;
		}
	}
	rep(i, 0, A + 1) dp[i][0] = 1;
	rep(j, 0, B + 1) dp[0][j] = 1;
	rep(i, 1, A + 1) {
		rep(j, 1, B + 1) {
			dp[i][j] = (dp[i - 1][j] + CS[j - 1][min(i, j - 1) + 1]) % mod;
		}
	}
	ll res = 0;
	rep(i, 0, A + 1) {
		rep(j, 0, B) {
			if(i + j > A) continue;
			res += dp[i][j] * C[B - 1][j];
			res %= mod;
		}
	}
	cout << res << "\n";
}

CODE FESTIVAL 2017 qual B,D: 101 to 010

http://code-festival-2017-qualb.contest.atcoder.jp/tasks/code_festival_2017_qualb_d

結局10111111...または111111...101となる文字列を見つける問題になるんですが、言い換えが甘くて場合分けをミスり、WAを量産しました。証明しながら解くことの重要性を感じた問題でした。

あとprevの情報を使うdpは覚えるものを複雑化するより、遷移を飛ばす感じにした方が実装楽そうですね。

int N;
string str;

int dp[MAX_N];
int nex[MAX_N];

void solve() {
	cin >> N >> str; str += "0";
	nex[N] = N;
	rer(i, N, 0) {
		if(str[i + 1] == '0') nex[i] = i + 1;
		else nex[i] = nex[i + 1];
		//debug(i, nex[i]);
	}
	rep(i, 0, N) {
		MAX(dp[i + 1], dp[i]);
		if(str[i] == '1') {
			if(nex[i] + 1 < N && str[nex[i] + 1] == '1') {
				MAX(dp[nex[i] + 2], dp[i] + nex[i] - i);
			}

			if(i + 1 < N && str[i + 1] == '0' && nex[i + 1] > i + 2) {
				MAX(dp[nex[i + 1]], dp[i] + nex[i + 1] - (i + 1) - 1);
				MAX(dp[nex[i + 1] - 1], dp[i] + nex[i + 1] - (i + 1) - 2);
			}
		}
	}
	cout << dp[N] << "\n";
}

CODE FESTIVAL 2017 qual A,E: Modern Painting

http://code-festival-2017-quala.contest.atcoder.jp/tasks/code_festival_2017_quala_e

まず最初に縦の向きに人を動かすとします。すると領域は二つに分断されます。そのうち一つについて横に進む人がX人、上から下に進む人がY人、下から上に進む人がZ人だとします。
するとその領域の塗り方はC(X+Y+Z-1,X-1)になります。なんでかは次の図を見てもらえればわかるかな…。

f:id:omochangram:20171007035342j:plain

このように一部反転させて、横に進む人が通過した領域を囲うようにして移動することを考えると、座標(0,0)から座標(X,Y+Z)までの、(Y座標)=Yの時には必ず一回は右に行くような移動方法の個数が領域の塗り方の個数と一致します。Y座標に初めてついたときのX座標をxとして[x,x+1)の区間を切り取れば、結局座標(X-1,Y+Z)へ移動する方法の個数を求めればいいことになり、C(X+Y+Z-1,X-1)となります。

これが求められたら後は累積和チックに求めれば計算量はO(N+M)となり十分高速です。

コンビネーションの複雑な式が出そうになったら、意味を考えて経路の問題に帰着させることを考えるといいかも?

上から下に進む人と右から左に進む人しかいなかったら簡単にできるなぁとは思いましたが、一回縦に分断した領域についてこんな簡単な式になるとは思いませんでした。

ll fac[MAX_N], inv[MAX_N], fiv[MAX_N]; //fiv:inv(fac(i))
ll pow2[MAX_N];

void C_init(int n) {
	fac[0] = fac[1] = 1; inv[1] = 1;
	fiv[0] = fiv[1] = 1;
	pow2[0] = 1; pow2[1] = 2;
	rep(i, 2, n + 1) {
		fac[i] = fac[i - 1] * i % mod;
		inv[i] = mod - inv[mod % i] * (mod / i) % mod;
		fiv[i] = fiv[i - 1] * inv[i] % mod;
		pow2[i] = pow2[i - 1] * 2 % mod;
	}
}

ll getC(int a, int b) { //assume a >= b
	if(a < b || a < 0 || b < 0) return 0;
	return fac[a] * fiv[b] % mod * fiv[a - b] % mod;
}

int N, M;
string A, B, C, D;
ll bitcnt[MAX_N];
ll sum[MAX_N];

ll solve2(int y, int z, int x) {
	if(x == 0) {
		if(y == 0 && z == 0) return 1;
		else return 0;
	}
	else return getC(x + y + z - 1, x - 1);
}

ll solve_sub() {
	int ccnt = accumulate(all(C), 0) - '0' * M;
	int dcnt = accumulate(all(D), 0) - '0' * M;
	int acnt = 0, bcnt = 0;
	bitcnt[0] = 0;
	rep(i, 0, N) bitcnt[i + 1] = bitcnt[i] + ((A[i] == '1' && B[i] == '1') ? 1 : 0);
	sum[N] = 0;
	rer(i, N, 0) {
		if(A[i] == '0' && B[i] == '0') sum[i] = sum[i + 1];
		else {
			sum[i] = sum[i + 1] + solve2(acnt, bcnt, dcnt) * pow2[bitcnt[i + 1]] % mod;
			sum[i] %= mod;
		}
		if(A[i] == '1') acnt++;
		if(B[i] == '1') bcnt++;
	}
	acnt = 0; bcnt = 0;
	acnt = 0; bcnt = 0;
	ll res = 0;
	ll mul = 1;
	rep(i, 0, N) {
		if(A[i] == '0' && B[i] == '0') continue;
		ll a = solve2(acnt, bcnt, ccnt);
		ADD(res, a * mul % mod * sum[i] % mod);
		if(A[i] == '1') acnt++;
		if(B[i] == '1') bcnt++;
		if(A[i] == '1' && B[i] == '1') MUL(mul, inv[2]);
	}
	return res;
}



void solve() {
	C_init(600010);
	cin >> N >> M;
	cin >> A >> B >> C >> D;
	int a = accumulate(all(A), 0) - '0' * N;
	int b = accumulate(all(B), 0) - '0' * N;
	int c = accumulate(all(C), 0) - '0' * M;
	int d = accumulate(all(D), 0) - '0' * M;
	if(a == 0 && b == 0 && c == 0 && d == 0) {
		cout << 1 << "\n";
		return;
	}
	ll ans = 0;
	if(a != 0 || b != 0) ADD(ans, solve_sub());
	if(c != 0 || d != 0) {
		swap(N, M);
		swap(A, C);
		swap(B, D);
		ADD(ans, solve_sub());
	}
	cout << ans << "\n";
}

100post目だやったぁ。