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 \ n) で求める。

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

(全体像)
公開鍵は、 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;
	}
}

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 -- )

			if (!(k & (1 << b))) dp[k|(1<<b)] = merge(dp[k|(1<<b)],dp[k]);
		}

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}) を求めるのは典型 ! 二分探索 + 貪欲

ARC094/ABC093 D:Worst Case

700点問題のD問題。むずかった。
D - Worst Case
高橋君のとった得点を A,B ( A \leq B) とする。
高橋君よりスコアの小さい最大の人数を求めたいので、1つ目のコンテストが高いスコアならば、2回目は低いスコア(逆もしかり)というように掛け算のペアを作ってあげる。(考察1)
(わからなければ解説放送へ)
www.youtube.com
高橋君よりスコアが小さい人のスコア a,b ( a \leq b) は必ず a < A または  b < B を満たすので
まず、 a < A の場合を考えてみる。
(考察1)より、以下のようにペアを考えることができる。
f:id:banboooo:20180408143002j:plain
すると、この場合はすべてのペアの積がABより小さくなることがわかる。
次に b < B の場合を考えてみる。先ほどのように 1 \sim (B-1) を元に、対応する A より大きい数を考えても、必ずABより大きいことが保証される数列が見つからない。
ここで発想を転換して、(A+1) \sim を元に、対応する数を考える。すると、積の最大値がABを超えない範囲で数列の長さを最大にする問題と捉えることができる。
ここまでを整理すると以下の図になる。
f:id:banboooo:20180408152026j:plain
ここからは2つほど解法があると思ったのだが、2つ目の方は実装で挫折してしまった。1つ目は、editorialの方法。2つ目は2分探索と3分探索を組み合わせる方法。
考え方などで間違っていたら教えていただきたいです。

1.editorialの方法
https://beta.atcoder.jp/contests/arc094/submissions/2325009
 C^2 < AB を満たす C を取ると、 C(C+1) または  C^2 が積の最大値となることを利用する方法。
これは掛け算のペアの和が常に一定ならば
ペアの数が同じくらいの大きさの時に、積の値が最大となるからである。
さらに例外的に、 A = B の時と  A+1 = B の時は、 C = A となってしまうので、別で処理しないといけない。
これを踏まえると以下のように考えることができる。


2.三分探索 & 二分探索の方法 (仮)
https://beta.atcoder.jp/contests/arc094/submissions/2325717 <- 不正解になってしまった提出。
 C の値をとることに気づかない場合は、こっち。
積の最大値を求めたいため、探索をするが逐次探索では当然間に合わない。
そこで、積の値の関数の形に注目してみると上凸の関数になっているはずである。

こういう時は、三分探索を使えば、極値、今回で言えば積の最大値を非常に効率よく求められる。
極値 Xの関数と見れば、この関数は単調非減少の関数なので、こちらは二分探索でXを定めてあげれば答えがもとまる。