AGC028 B問題 Removing Blocks

beta.atcoder.jp

editorialで確率がどうのこうの書いてあって ??? となったので。

方針

最終的に求めたい答えは当たり前だけど  a_0 A_0 + a_1 A_1 + \cdots + a_n A_n の形でかける.
そして ,  a_i の値を求めにいく.

  •   a_i = \displaystyle \sum_{j} ( A_j を取り除く時にA_iが連結である場合の数 )  で求めること.
  •  A_j を取り除く時にA_iが連結である場合の数を全体の  N! 通りの中の割合で求めること.

この二つが大きなポイントだと思う. 二つ目の方から考える.

2つ目のポイントの考察

※下記では j \leq i として議論する.
 A_j を取り除く時にA_iが連結である  \Leftrightarrow  A_j , \cdots , A_i 間で A_jが最初に除かれる.
と言い換えることができる.では,その場合の数はN!のうちどのくらいの割合なのか?

 A_j \cdots A_i 間で,初めて取り除かれるものを A_kとし , これによって事象を定義する.
すなわち , 全事象  \Omega = \displaystyle \bigcup_{k=j}^i   \left\{ 事象k ; A_j \cdots A_i 間でA_kが初めて除かれる \right\}
これらは排反の事象であり , 確率(場合の数の割合)はどの事象も等しく  P(k) = \dfrac{1}{i - j + 1} となっている.

A_jが最初に除かれる場合の数を知りたかったので, 全体の場合の数をかけて

 N! P(j) = \dfrac{N!}{i - j + 1}

あと ,  j \leq i と条件を取っ払えば

 N! P(j) = \dfrac{N!}{|i - j| + 1}

1つ目のポイントの考察

あとはさっきの結果を代入してあげればok

 a_i = \displaystyle \sum_{j}  \dfrac{N!}{|i - j|+ 1}

最後にもう一工夫.
  0 \leq |i-j| \leq N より ,  \displaystyle \sum_{k = 0}^N \dfrac{N!}{k+ 1} を先に計算(累積和)しておくことができる.  O(N)
すると  a_i  O(1)で求まるようになる.

なので,全体は答え  a_0 A_i + a_1 A_1 + \cdots + a_n A_n O(N)で求まる.

ソースコード

N = int(input())
A = list(map(int,input().split(" ")))
mod = 10**9+7

# N!計算
fact = 1
for i in range(2,N+1):
    fact *= i
    fact %= mod

# 場合の数計算
rev = [ (pow((i),mod-2,mod)*fact)%mod for i in range(1,N+1) ]
# 累積和
for i in range(N-1):
    rev[i+1] += rev[i]
    rev[i+1] %= mod

# a_0 A_0 + a_1 A_1 + ... + a_n A_n を計算
ans = 0
for i in range(N):
    ans += A[i] * (rev[i] + rev[N-1-i] - rev[0]) % mod
    ans %= mod

print(ans)

感想

解けてから見ると問題文にわざわざ N!通りとか書いてくれてるのはヒントだったのかなぁと思ったり。

ARC102 D問題 All Your Paths are Different Lengths

なんか最近難化傾向にあるような...
D - All Your Paths are Different Lengths

<考察>
1. 以下のようにノード  n , m 間に長さ  0 , 2^{n} の二つの辺を貼ると ,
  L = 2 ^ {node} の時のグラフが表現できる.

f:id:banboooo:20180902143647j:plain:w300

2.  [ 0 , L ) をサイズが2の累乗の区間の和で表すことを考える.
  ( ルール : サイズの大きい区間から貪欲に決めていくこと )

  例1 :
    [0 , 10) =  [0 , 8) \  \bigcup \ [8 , 10) \\ 
   = [0 , 8 )\  \bigcup \  [ 8 + 0 , 8 + 2 )
  例2 :
    [0 . 15) =  [0 , 8)\  \bigcup \  [8 , 12) \  \bigcup \ [12 , 14) \ \bigcup \ [14 , 15) \\ 
                                = [0 , 8) \ \bigcup \ [ 8 + 0 , 8 + 4 ) \ \bigcup \ [ 12 + 0 , 12 + 2 ) \ \bigcup \ [ 14 + 0 , 14 + 1 )

  これと同じようにグラフをそれぞれ作ると ...
f:id:banboooo:20180902152739j:plain:w500

こうすると , 同じ長さのパスを一つも作ることなく , グラフを構築することができました .
pythonでの実装

L = int(input())
# ノード数を決める
for i in range(1,21):
    if (1 << i) > L:
        node = i
        break

# 必ず 0と2^Nの辺を貼る
edge = []
for i in range(node-1):
    edge.append((i+1,i+2,0))
    edge.append((i+1,i+2,1<<i))

# 2の累乗の区間に分割
for i in range(node-1):
    if (L >> i) & 1:
        edge.append((i+1,node,L-(1<<i)))
        L = L-(1<<i)

#出力
print(str(dig) + " " + str(len(edge)))
s = "{0} {1} {2}"
for v,v2,w in edge:
    print(s.format(v,v2,w))


  
  



RSA暗号のpython実装

競プロとは関係ないですが。授業で聞いて面白かったので実装してみました。
アルゴリズムの解説

(準備)

  1. 素数  p1,p2 を決める (十分大きな数)
  2.  L = lcm(p1-1,p2-1) を求め、L と互いに素な数  e を1つランダムに選ぶ。
  3.  e x \equiv 1 (mod L ) となる  x ( x > 0) を求める。(  ex + Ly = 1 の特殊解を求めることと同値 )

(暗号化)
平文  S (数字) に対し、暗号  C  C = S^e (mod \ p1\times p2) で求める。

(復号化)
暗号  C (数字) に対し、元の平文  S C^{x} (mod \ p1\times p2) で求める。

(全体像)
公開鍵は、 p1 \times p2  e である。
 p1 , p2 が十分に大きい時、 p1 \times p2 から  p1 , p2を求めることが現実的でないことを利用している。

from random import randint
class RSA:
	def __init__(self,p1,p2):
		self.__p1,self.__p2,self.__d,self.e,self.n = self.__makeKey(p1,p2)

	def __gcd (self,a,b):
		if a < b:
			a,b = b,a
		return a if b == 0 else self.__gcd(b,a%b)

	def __lcm (self,a,b):
		return a*b//self.__gcd(a,b)

	def __extended_euclid(self,a,b,k = 1):
		if k != 1:
			g = self.__gcd(a,b)
			a //= g
			b //= g

		cc = gg = 1
		ee = ff = 0
		while b:
			div,mod = divmod(a,b)
			hh = cc - div * ee
			ii = ff - div * gg
			a,b = b,mod
			cc,ee = ee,hh
			ff,gg = gg,ii

		return (cc,ff)

	def __choose_e(self,n,e = 0):
		if e != 0:
			return e
		while True:
			k = randint(n+100,n+1000)
			if self.__gcd(k,n) == 1:
				return k

	def __makeKey (self,p1,p2):
		n = p1*p2
		l = self.__lcm(p1-1,p2-1)
		e = self.__choose_e(l)
		(d,se) = self.__extended_euclid(e,l)
		if d <= 0:
			d += l 
		return p1,p2,d,e,n

        # 暗号化(数字)
	def encode (self,x):
		if x > self.n:
			print("x must be under {0}".forat(self.n))
			exit(1)
		return pow(x,self.e,self.n)
       # 復号化(数字)
	def decode (self,c):
		return pow(c,self.__d,self.n)
       # 暗号化(文字)
	def encode_str (self,s):
		ret = []
		for i in range(len(s)):
			ret.append(self.encode(ord(s[i])))
		return ret
        # 復号化(文字)
	def decode_str (self,s):
		ret = []
		for i in range(len(s)):
			ret.append(chr(self.decode(s[i])))
		return "".join(ret)

実行例

A = RSA(149,151) # 149 と 151

# 公開鍵
print (A.e,A.n)

# 数字の暗号化
Code = A.encode(57)
print(Code)

#数字の復号化
Raw = A.decode(Code)
print(Raw) # 57に戻る

#文字の暗号化
Code = A.encode_str("Hello")
print(Code)

#文字の復号化
Raw = A.decode_str(Code)
print(Raw) # "Hello"に戻る

ARC100 E問題 OrPlusMax (高速ゼータ変換)

E - Or Plus Max
高速ゼータ変換の問題。といっても貼るだけの問題ではなく、集合や演算を自分で考察してうまくテンプレに落としてやる必要がある。
また、やってることはbitDPなので言葉を知らなくても解ける(僕はコンテスト中には解けてない)
この記事では、どこの部分がいわゆる高速ゼータ変換のテンプレになっているかを書いておく。

  1. そもそも高速ゼータ変換とは ?

こちらの記事を参考にさせていただきました。
https://topcoder.g.hatena.ne.jp/iwiwi/20120422/1335065228


高速ゼータ変換
集合に対する関数  f(S) に対し、全集合  S に対して、 g(S) = \sum_{S \subseteq T} f(T)

これを実装したものが

 rep (i, n) rep (b, 1 << n) if ((b & (1 << i)) == 0 )  a[b] += a[b | (1 << i)];

と書けるそう。

  1. 本問ではどのように使われているのか ?

結論から先に言うと、

  • 集合  S := \{ i | 数nのiビット目が0 \}
  • 演算  merge := 1番目に大きい数と2番目に大きい数を更新

と定義する。ソースコードはこんな感じ。

 rep (i, n) rep (b, 1 << n) if ((b & (1 << i)) == 1 )  
        a[b] = merge( a[b] , a[b ^ (1 << i)] );

( おまけ )
集合を S := \{ i | 数nのiビット目が1 \} と取ることもできるが、
この時は  g(S) = \sum_{S \supseteq T} f(T) になっており、
理解の足りない僕はあれ?包含関係逆じゃね ? いいの ? とか考えてしまった。
また、こちらの方針で実装する場合は以下のようになる。

for (int i = (n-1); i >= 0 ; i-- )
    for (int b = (1 << n) - 1 ; b >=  0 ; b-- )
        if (!(b & (1 << i))) dp[b|(1<<i)] = merge(dp[b|(1<<i)],dp[b]);

ループの向きを逆にしないといけない。これは大いにバグの原因になりそう。

最後に僕のE問題のソースコード

//演算
#include <bits/stdc++.h>
vector<long> merge(vector<long> &A,vector<long> &B) {
	vector<long> cp(A);
	REP(i,B.size()) cp.push_back(B[i]);
	sort((cp).begin(),(cp).end(),greater<long>());
	return vector<long>(cp.begin(),cp.begin()+min(2,(int)cp.size()));
}

int main(void) {
	int N;
	cin >> N;
	vector<long> A(1<<N);
	for(long i = 0;i < (1<<N);i++) cin >> A[i];
	vector<vector<long>> dp(1<<N,vector<long>(2,0));
	for(long i = 0;i < (1<<N);i++) dp[i] = {A[i]};

	/*ゼータ変換部分*/
	for(int i = 0;i < N;i++) REP(b,1<<N) if ((b & (1 << i))) {
		dp[b] = merge(dp[b],dp[b ^ (1 << i)]);
	}

	long ans = 0;
	for(int i = 1;i < (1<<N);i++) {
		ans = max(ans,dp[i][0]+dp[i][1]);
		cout << ans << endl;
	}
}

codeFlyer C問題 徒歩圏内

解けなくて悔しいので記事に。考察の流れとかも書いておく。
排反をとるというのがtwitterで見ててわかりやすそうだったので実装。

bitflyer2018-qual.contest.atcoder.jp


考察した(ダメ)解法
1.  O(N^2 \log N ) : 左端と真ん中を決めて右端の数を二分探索で数える。左端と真ん中の決め方をしゃくとりできるかなと思ったが無理。{}_n C _2通り試さないとだめ。当然TLE。

2. 左端と右端を決めて、真ん中を二分探索で数える。1と本質的には変わらず。

二分探索と累積和を使う O(N \log N) 解法
正しい解法
左端の都市、真ん中の都市、右端の都市のindexをそれぞれ  i , j , k ( i < j < k ) と書く。
全ての  ( i , j , k ) の組み合わせは、 {}_n C _3 通り
ここから条件に合わないものを引くことを考える。

左端の都市  i を順に見ていく。

1.  (X_j - X_i) > D の場合 ->  {}_(X_iから右にDより離れている都市数) C _2 通り
真ん中が左端からDより離れている場合、右端も左端からDより離れている。
よって、(左端から右にDより離れている都市) の中から2つ選ぶ組み合わせに対応する。

2.  (X_ j - X_i ) \leqq D かつ ( X_k - X_j ) > D の場合 ->  \sum_{j}{} ( X_j から右にDより離れている都市数)
真ん中は左端からD以内だが、右端が真ん中から D以上離れている。
真ん中を決めうちし、右端の候補数を足し合わせたい。

3.  (X_k - X_i ) \leqq D の場合 ->  {}_(X_iから右にD以内にある都市数) C _2 通り
左端と右端がD以内の距離にある場合。この時、真ん中の点もD以内にある。
よって、(左端から右にD以内の都市) の中から2つ選ぶ組み合わせに対応する。

1,2,3で全ての場合を尽くしているので、これを全ての i に対して合計して、最後全体から引けば良い。

ここで、1,2で何度も  ( Dより離れている都市数 ) を使っており、2ではこれの累積値を取っていることから、
あらかじめ累積和を取っておくと2がスムーズに計算できそう。
各都市に対して、右にD以上離れている最小のindexの都市を二分探索で求める。
-> 累積和をとる。累積前も取っておくといいかな (?)

3の  (D以内の都市数) (Dより離れている都市数) から簡単に求まるので心配なし。

ソースコード
Submission #2603957 - codeFlyer (bitFlyer Programming Contest) | AtCoder

自作ライブラリ(行列編)

気が向いたので作ってみました。固有値求める機能とかも追加したいけどとりあえず後回し。

四則演算,累乗,ビット演算,転置,行列式,逆行列とかその辺は実装しました〜

ABC034 ~食塩水~

D - 食塩水
解いたのでメモ。悩んだけど、どうやら典型の模様。
蟻本p132の「平均最大化」とほとんど同じ問題。
最初は、ナップザック問題っぽいなぁと思って、DPを必死に考えてたけど
よくよく考えてみれば以下の式は別物。
本問 :  max(\dfrac{\sum v}{\sum w})
ナップザック :  max(\sum v) \ \ {\sum w \leq const}

~ アルゴリズム ~
目標の  p を決めてあげて、この濃度以上のものが作れるか判定する。
濃度  p 以上のものが作れる or 作れないの境界値が今回の答えの値。

まず、判定をどうするか。目標を  P^* とすれば各容器  i に対し、 w_i \times (P^* - p_i) を求めて、貪欲に大きいものからとって 0以上なら  P^* は作れる。0未満なら作れないとやれば良い。
 (実際の食塩の量) = p_i w_i + p_{i+1} w_{i+1} + ...
 (必要な食塩の量) = P^* \sum w = P^*w_i + P^*w_{i+1} +...

判定ができれば、あとは境界値を探れば良い。

ソースコード
Submission #2374269 - AtCoder Beginner Contest 034


まとめ :
 max(\dfrac{\sum v}{\sum w}) を求めるのは典型 ! 二分探索 + 貪欲