Suffix Arrayを使ってみた
全く解き方が分からなかったMountain HolidaysのEditorialが出たので解いてみました．
Mountain Holidays
整数の配列height(1 <= |height| <= 100000)が与えられる．このとき，連続部分列の数を数えよ．
次の条件のいずれかを満たすとき(i1, j1)と(i2, j2)の2つの連続部分列は異なる．
1. j1 – i1 != j2 – i2
2. height[i1+k] – height[i1+k-1] != height[i2+k] – height[i2+k–1] を満たすkが (1, 2, ... j1-i1)の範囲に存在．
この問題は異なる部分文字列の個数を数える問題に帰着できます．height[i] - height[i-1]をその整数に対応する文字と考えます．
しかし，この問題では|height|が大きいため，ハッシュなどを使用して高速な探索を行ったとしても部分文字列の個数を探すのにO(N^2)かかってしまいTLEします．
そこで，heightのsuffixを辞書順に並べるSuffix Arrayを構築します．Suffix Arrayを用いることでi番目のsuffixとi+1番目のsuffixのLCP(Longest Common Prefix)を高速に求めることができ，（全ての部分文字列の個数） - （i番目のsuffixとi+1番目のsuffixのLCPの総和）で異なる部分文字列の個数を出すことが出来ます．
Suffix Arrayの構築は普通にやると，ソートにO(NlogN)，1回の比較にO(N)かかるため，O(N^2logN)となり間に合いません．
そこで，2^iの長さごとに辞書順を求めておき，2^(i+1)の長さの辞書順を求める際に2^iの長さの辞書順を用いることでO(N(logN)^2)でSuffix Arrayを得ることが出来ます．Suffix Arrayを得ることが出来れば1つのLCPをO(logN)で求められるため，LCPの総和をO(NlogN)で求めることが出来，TLEせずに実行することが出来ます．
#include <iostream> #include <cstdio> #include <cstdlib> #include <cmath> #include <vector> #include <algorithm> #include <string> #include <map> #include <set> #include <queue> #include <stack> #include <climits> #include <sstream> #include <functional> #include <complex> using namespace std; #define len(array) (sizeof (array) / sizeof *(array)) #define rep(i, s, e) for(int i = s;i < e;i++) #define Rep(i, e) for(int i = 0;i < e;i++) #define rrep(i, e, s) for(int i = e;s <= i;i--) #define Rrep(i, e) for(int i = e;0 <= i;i--) #define mrep(i, e, t1, t2) for(map<t1, t2>::iterator i = e.begin(); i != e.end(); i++) #define vrange(v) v.begin(), v.end() #define vrrange(v) v.rbegin(), v.rend() #define vsort(v) sort(vrange(v)) #define vrsort(v) sort(vrrange(v)) #define arange(a) a, a + len(a) #define asort(a) sort(arange(a)) #define arsort(a, t) sort(arange(a), greater<t>()) #define afill(a, v) fill(arange(a), v) #define afill2(a, v, t) fill((t *)a, (t *)(a + len(a)), v) #define fmax(a, b) (a < b? b : a) #define fmin(a, b) (a > b? b : a) #define fabs(a) (a < 0? -a : a) #define pb push_back #define rg(i, s, t) s <= i && i < t //#define X real() //#define Y imag() //typedef unsigned int ui; typedef long long ll; typedef unsigned long long ull; typedef pair<int, int> P; //typedef complex<double> p; const int INF = (int)2e9; const int MOD = (int)1e9 + 9; const double EPS = 1e-10; //const int dx[] = {1, -1, 0, 0, 1, -1, -1, 1}; //const int dy[] = {0, 0, 1, -1, -1, -1, 1, 1}; //const int weight[] = {0,1,10,100,1000,10000,100000,1000000,10000000}; //priority_queue< int, vector<int>, greater<int> > q; //Suffix Array //http://www.stanford.edu/class/cs97si/suffix-array.pdf //http://discuss.codechef.com/questions/18105/mou1h-editorial #define MAX_N 100005 #define MAXLG 18 #define MIN_VALUE -INF struct entry { int nr[2], p; bool operator <(const struct entry &e) const{ return nr[0] == e.nr[0] ? (nr[1] == e.nr[1]? p < e.p : nr[1] < e.nr[1]) : (nr[0] < e.nr[0]); } bool operator >(const struct entry &e) const{ return nr[0] == e.nr[0] ? (nr[1] == e.nr[1]? p > e.p : nr[1] > e.nr[1]) : (nr[0] > e.nr[0]); } } L[MAX_N]; int SA[MAXLG][MAX_N], N, suffix[MAX_N], LCP[MAX_N]; int A[MAX_N], D[MAX_N]; //SA[k][i] : the order of the subsequence of A of length 2^k starting at position i. //return the last row index of matrix SA // O(N (logN)^2) int makeSuffixArray(int n) { // cin >> A; int stp = 1, cnt = 1; N = n; Rep(i, N) SA[0][i] = D[i]; while(cnt < N){ Rep(i, N){ L[i].nr[0] = SA[stp - 1][i]; L[i].nr[1] = (i + cnt < N)? SA[stp - 1][i + cnt] : MIN_VALUE; L[i].p = i; } sort(L, L + N); Rep(i, N){ SA[stp][L[i].p] = (i > 0 && L[i].nr[0] == L[i - 1].nr[0] && L[i].nr[1] == L[i - 1].nr[1])? SA[stp][L[i - 1].p] : i; } stp++; cnt *= 2; } return stp-1; } //Computing the longest common prefix (LCP) // O(log N) int lcp(int x, int y, int stp) { int ret = 0; if (x == y) return N - x; for (int k = stp; k >= 0 && x < N && y < N; k--){ if (SA[k][x] == SA[k][y]){ x += (1 << k); y += (1 << k); ret += (1 << k); } } return ret; } //結果はsuffixに格納． int getSuffixArray(int n){ int idx = makeSuffixArray(n); Rep(i, N){ suffix[SA[idx][i]] = i; } return idx; } //結果はLCPに格納 void getLCPArray(int stp){ LCP[0] = 0; rep(i, 1, N){ LCP[i] = lcp(suffix[i-1], suffix[i], stp); } } void doIt(){ int t, n; ll ans, sum; scanf("%d", &t); while(t--){ scanf("%d", &n); Rep(i, n) scanf("%d", &A[i]); if(n == 1){ printf("0\n"); continue; } else if(n == 2){ printf("1\n"); continue; } ans = n; ans *= n-1; ans /= 2; ans %= MOD; n--; rep(i, 0, n) D[i] = A[i+1] - A[i]; int idx = getSuffixArray(n); getLCPArray(idx); // printf("idx = %d\n", idx); // printf("idx2 = %f\n", ceil(log(N)/log(2))); Rep(i, n){ ans -= LCP[i]; // printf("LCP[%d] = %d\n", i, LCP[i]); if(ans < 0) ans += MOD; ans %= MOD; } printf("%lld\n", ans); } } int main() { doIt(); return 0; }
この問題，全く通せなくてすごい苦労しました．Editorial通りに実装しているし，他の人のコードを見ても同じようにやっているのになんで通らないんだ？と数時間格闘し，最後は自力は諦めて通っている人のコードと差分を取りながらデバッグ．
その結果，MODに使う数字がずれていたと判明（！！？）．多くの問題ではMODは10^9+7なのですが，この問題は10^9+9というしょうもないオチでした．これはひどい．脱力・・・．
参考
http://discuss.codechef.com/questions/18105/mou1h-editorial
http://www.stanford.edu/class/cs97si/suffix-array.pdf
SRM609参加記
初めてdiv2の本番で全完出来てテンションが上がったので書いちゃいます．
MagicalStringDiv2(250)
'>'と'<'からなる偶数長の文字列が与えられる．これを'>''>'...'<''<'の形にするには最低何文字変換しなければならないか．
単に数えるだけ．
class MagicalStringDiv2 { public: int minimalMoves(string S) { int len = S.length(); int ret = 0; Rep(i, len){ if(i < len / 2){ if(S[i] == '<') ret++; } else{ if(S[i] == '>') ret++; } } return ret; }
PackingBallsDiv2(500)
赤，緑，青のボールがそれぞれR，G，B個ある．それらを使って可能な限り少ないパッケージを作る．パッケージは
1. 1〜3個からなる
2. 全て同じ色か，全て異なる色からなる
の条件を満たさなければならない．
なんか前にこどふぉで解いた問題に似てるなあと思いつつ考察．
こどふぉの時は変な条件を考えて条件分岐したら死んだので答え列挙系にしたい．
そこで，試しに全部列挙してみたら最大ケースでもいけそうだったので提出．
しかし，ループ毎に変数の値を戻し忘れるという変なミスをしたせいで再提出．
class PackingBallsDiv2 { public: int minPacks(int R, int G, int B) { int ans = INF; Rep(i, 101){ int res = 0; int rr = R, gg = G, bb = B; if(R >= i && G >= i && B >= i){ res += i; rr -= i; gg -= i; bb -= i; Rep(j, 101){ if(rr >= j && gg >= j){ // printf("%d, %d\n", i, j); res += j; rr -= j; gg -= j; Rep(k, 101){ if(rr >= k && bb >= k){ res += k; rr -= k; bb -= k; Rep(l, 101){ // printf("%d, %d, %d, %d, gg = %d, bb = %d\n", i, j, k, l, gg, bb); if(gg >= l && bb >= l){ // printf("%d, %d, %d, %d\n", i, j, k, l); res += l; gg -= l; bb -= l; int rrr = rr, ggg = gg, bbb = bb, rres = res; rrep(m, 3, 1){ rres += rrr / m; rres += ggg / m; rres += bbb / m; rrr %= m; ggg %= m; bbb %= m; } ans = min(ans, rres); res -= l; gg += l; bb += l; } } res -= k; rr += k; bb += k; } } res -= j; rr += j; gg += j; } } res -= i; rr += i; gg += i; bb += i; } } return ans; }
VocaloidsAndSongs(950)
ボーカロイドのGumi，Ia，MayuはS曲からなるアルバムを作ることにした．
曲は1人〜3人が歌う．全ての組み合わせ数をXとしたとき，X mod 1,000,000,007を求めよ．
残り30分くらいだったのでさすがに厳しいかなーとオープンしたHardがボカロ題材だったのでちょっと嬉しかった（@semiexpさんの問題だったらしい．これからもボカロの問題を出してくれるといいなあ・・・）．しかも問題文短い．
はじめ，「全部列挙して歌ってる人数が0の曲が含まれている場合を引けばいけそう！」と思い，包除原理とか使うのかなーと考えてましたが，「あれ？これってDPじゃ無理なのか？」と思い至り，試しに実装してみたらテストケースパスしたので提出．
dp[i曲め][GUMIの歌った数][IAの歌った数][MAYUの歌った数]でDPします．
ll dp[52][52][52][52]; class VocaloidsAndSongs { public: int count(int S, int gumi, int ia, int mayu) { afill2(dp, 0, ll); dp[0][0][0][0] = 1; rep(i, 1, S+1){ Rep(j, gumi+1){ Rep(k, ia+1){ Rep(l, mayu+1){ if(i && j) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j-1][k][l]) % MOD; if(i && k) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j][k-1][l]) % MOD; if(i && l) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j][k][l-1]) % MOD; if(i-1>=0 && l && j) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j-1][k][l-1]) % MOD; if(i-1>=0 && l && k) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j][k-1][l-1]) % MOD; if(i-1>=0 && j && k) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j-1][k-1][l]) % MOD; if(i-1 >= 0 && l && j && k) dp[i][j][k][l] = (dp[i][j][k][l] + dp[i-1][j-1][k-1][l-1]) % MOD; } } } } return dp[S][gumi][ia][mayu]; }
初めて全完でチャレンジフェイズを迎えましたが，SRMは自信があってもよく落ちるので「どうせ落ちるんだろうなー」と思ってました．
チャレンジフェイズではMedで落とせそうなコードを発見．数行で書かれてるし，なんかおかしそう．でもPythonで書かれていたのでよく分からない．
そこで写経．しかし，初めてemacsでPythonを書きましたがオートインデントは効かないし，BackSpaceによるインデントの削除は出来ないし四苦八苦（僕のemacsだけ？）．
ようやく書けたものの結局落とせるケースはよく分からず（このコードはSystem Testで落ちてた）．
結局244.30，191.00，646.06の1081.36でRoom2位，全体72位でした（初めての全完で高まっていたけどそこまでいいわけでもなかった）．
1074 -> 1146
もしかしてDiv1戻れるかなーとか期待していましたが普通に無理でした．
Div1への道は本当に険しい・・・．
ARC017に参加しました
ようやく修論発表も無事（？？？）終わり、少し余裕が出来たのでARC017に参加。今回こそはCを解くと意気込んでいました。
A - 素数、コンテスト、素数
N(17≤N≤1,000,000)が与えられるので、Nが素数かどうか判定する。
コピペをやるだけ。
始めsqrt(1,000,000)とか考えてたけど、冷静に考えてO(n)で普通にいけた。
B - 解像度が低い。
要素数 N(1≤N≤300,000)の配列Aと整数K(1≤K≤N) が与えられる。Aの中でK個連続で単調増加になっている部分は何個あるか？
A[i]番目までに単調増加になっている個数seq[i]を覚えておく。A[i]番目とA[i+1]番目を見て、A[i+1]の方が大きかったらその数に1を足し、そうでなかったら1をseq[i+1]とする。あとはK<=seq[i]となる個数を数えれば良い。
ここまでは簡単なのですが、僕の実装方法が悪くてn=1とk=1のコーナーケースが発生しており、2WA。
C - 無駄なものが嫌いな人
大きさ X(1≤X≤109) のナップサックと N(1≤N≤32) 個の品物がある。i 番目の品物の大きさは wi(1≤wi≤5*10^7) である。N 個の品物から、大きさの総和が無駄なくぴったり X となる選び方が何通りあるか？
始め、再帰で枝狩りしまくればいけるかと思ってやったけどさすがに2^32は強く普通にTLE。
終了前10分くらいで2^16ならいけるんじゃ？と半分全列挙を思いつく。
しかし、実装出来る訳もなく終了。
意気込んでいたCは解けず、88位と微妙な順位で終わってしまいました。
コンテスト終了後にCの半分全列挙を実装してみました。
半分全列挙とは全パターンを列挙することができないようなサイズの問題に対し、半分ずつ列挙することで全ての列挙を可能にするアルゴリズムです（from 蟻本）。
今回の場合、2^32は無理ですが、半分ずつ列挙すると2^16 * 2で半分からなる組み合わせが全て列挙でき、あとは片方の要素に対し二分探索を使うことで時間内におさめることが出来ます。
#include <iostream> #include <cstdio> #include <cstdlib> #include <cmath> #include <vector> #include <algorithm> #include <string> #include <map> #include <set> #include <queue> #include <stack> #include <climits> #include <sstream> #include <functional> #include <complex> using namespace std; #define len(array) (sizeof (array) / sizeof *(array)) #define rep(i, s, e) for(int i = s;i < e;i++) #define Rep(i, e) for(int i = 0;i < e;i++) #define rrep(i, e, s) for(int i = e;s <= i;i--) #define Rrep(i, e) for(int i = e;0 <= i;i--) #define vrange(v) v.begin(), v.end() #define vrrange(v) v.rbegin(), v.rend() #define vsort(v) sort(vrange(v)) #define vrsort(v) sort(vrrange(v)) #define arange(a) a, a + len(a) #define asort(a) sort(arange(a)) #define arsort(a, t) sort(arange(a), greater<t>()) #define afill(a, v) fill(arange(a), v) #define afill2(a, v, t) fill((t *)a, (t *)(a + len(a)), v) #define fmax(a, b) (a < b? b : a) #define fmin(a, b) (a > b? b : a) #define fabs(a) (a < 0? -a : a) #define pb push_back #define rg(i, s, t) s <= i && i < t //#define X real() //#define Y imag() //typedef unsigned int ui; typedef long long ll; typedef unsigned long long ull; typedef pair<int, int> P; //typedef complex<double> p; const int INF = (int)2e9; const int MOD = (int)1e9 + 7; const double EPS = 1e-10; //const int dx[] = {1, -1, 0, 0, 1, -1, -1, 1}; //const int dy[] = {0, 0, 1, -1, -1, -1, 1, 1}; //const int weight[] = {0,1,10,100,1000,10000,100000,1000000,10000000}; //priority_queue< int, vector<int>, greater<int> > q; // #define MAX_N 1000 void doIt(){ int n, x, w[33]; cin >> n >> x; Rep(i, n){ cin >> w[i]; } // sort(w, w+n, greater<int>()); map<int, int> grpA[2], grpB[2]; int prev = 0, now = 1, ptr; grpA[prev][0]++; map<int, int>::iterator it; Rep(i, n/2+1){ it = grpA[prev].begin(); grpA[now].clear(); while(it != grpA[prev].end()){ grpA[now][(*it).first + w[i]] += (*it).second; grpA[now][(*it).first] += (*it).second; it++; } prev = 1 - prev; now = 1 - now; } ptr = prev; prev = 0; now = 1; grpB[prev][0]++; rep(i, n/2+1, n){ it = grpB[prev].begin(); grpB[now].clear(); while(it != grpB[prev].end()){ grpB[now][(*it).first + w[i]] += (*it).second; grpB[now][(*it).first] += (*it).second; it++; } prev = 1 - prev; now = 1 - now; } it = grpA[ptr].begin(); ll ans = 0; while(it != grpA[ptr].end()){ // printf("itA = %d\n", (*it).first); if(grpB[prev].count(x - (*it).first)) ans += (*it).second * grpB[prev][x - (*it).first]; it++; } cout << ans << endl; } int main() { doIt(); return 0; }
上位の方々のコードを参考に実装。
今まで知らなかったんですが、map
CodeChef October Challenge 2013を頑張ってみた
October Challenge 2013に参加しました。
初めて7問（正確には6.143問）解けたので参加記書きます。以下は解いた人が多い順です。
Maxim and Dividers
10進数表記で4か7を含むものをoverluckyと呼ぶ。n(1 <= n <= 10^9)が与えられたとき、nを割り切れるoverluckyはいくつあるか。
sqrt(n)まででnを割り切れ、かつ4か7を含むものを数える。
Little Elephant and Bamboo
長さnの配列HとDが与えられる。配列Hに対し、1つの要素を1減らし、後の要素を1増やす操作を行うことが出来る。このとき、HをDと同じ配列にするのに何回操作を行えばよいか出力せよ。無理なら-1を出力せよ。
考え方を変えると、上記の操作は「全ての要素に1を足し、1つの要素から2を引く」と考えることが出来る。なので、もしHをDに変形可能ならば（Dの総和）-（Hの総和）が(n-2)で割り切れるのが必要条件。あとは個々の要素を見ていく。
Yet Another Nice Girl
長さnの正数の配列が少年と少女に与えられる。少年と少女はランダムに配列から要素を取り出して、勝負する。少年の値をx、少女の値をyとしたとき、x^y > y^xなら少年は少女にキスしてもらえる。少年がキスしてもらえる回数の期待値を出力せよ。
xを固定したとき、yが何ならx^y > y^xになるかを考えてみると、
(i) x = 1
yが何でも勝てない。
(ii) x = 2
y = 2, 3以外は勝てる。
(iii) x = 3
3以外勝てる。
(iv) x = その他
自分よりも大きい数字と1に勝てる。
なので、少女の配列をソートして、二分探索をかければ固定したxが勝てる期待値を求めることが出来る。少年の配列の全ての値に対して期待値を求めればよい。
Sereja and Transformation
長さnの配列aが与えられる。b[i] = min(a) - a[i] + sum(a)の変形をk回行った配列をrとする。q(r) = max(r)-min(r)とする。1 <=a[i] <= mであるとき、a[i]は何通り存在するか出力せよ(mod 10^9+7)。
変形を見て難しそうに思えるが、実は配列内の要素の相対的な差は変わらない。なので、配列aで最大値と最小値を固定し、その数を数えればよい。
Edit Distance on Grid
N*Mの二次元格子が与えられる。各格子は白か黒で塗られている。
(i) 辺を共有している格子は1秒で取り替えることができる。
(ii) 白い格子を黒くするのはC2秒で出来る。
(ii) 黒い格子を白くするのはC3秒で出来る。
全て白い格子にするか、全ての黒の格子が繋がっている格子にするには最短何秒必要か。
CodeChefに毎回1問くらいあるマラソン系問題。マラソンめっちゃ苦手なので、とりあえず白と黒の少ない方を反転させるコードを提出。
Kamehameha
二次元平面に魔物のいる座標がN個与えられる。悟空はかめはめ波で一行、一列の魔物を全て消し去ることが出来る。最小何発のかめはめ波で魔物を全て消し去れるか。
全く分からなかった問題。あまりに思いつかないので「まさか枝狩りしまくればいけるのか・・・？」と2^1000に挑むも当然TLE。
なんと蟻本にほぼ同じ問題が載っていました。最小辺被覆に帰着できるらしいです。二部グラフにおいては最小辺被覆の辺数 = グラフの最大マッチングが成り立つので、最大マッチングを求めればOKです。
Rating 3892 -> 4517
今回は蟻本のおかげでレートが大きく上がりました。しかし、CodeChefのレーティングシステムはよく分からない。
[図形][C++] 多角形クラスを作ってみた
多角形の面積を簡単に出す方法を@hirokazu1020さんに教わったので多角形クラスPolygonを実装してみました。ほとんど蟻本のコピーです。
//二次元座標(from 蟻本) typedef struct _PT{ double x, y; _PT() {} _PT(double x, double y) : x(x), y(y){ } _PT operator + (_PT p){ return _PT(x + p.x, y + p.y); } _PT operator - (_PT p){ return _PT(x - p.x, y - p.y); } _PT operator * (double d){ return _PT(d*x, d*y); } double dot(_PT p){ //pとの内積 return x * p.x + y * p.y; } double det(_PT p){ //pとの外積 return x * p.y - p.x * y; } double dist(_PT p){ //pとの距離の2乗 return (x-p.x)*(x-p.x) + (y-p.y)*(y-p.y); } bool operator <(const struct _PT &e) const{ return x == e.x? (y < e.y) : x < e.x; } bool operator >(const struct _PT &e) const{ return x == e.x? (y > e.y) : x > e.x; } class Polygon{ private: public: vector<pt> pts; int n; Polygon(){} Polygon(vector<pt> v){ n = v.size(); pts = v; } //凸包（蟻本ver） vector<pt> convex_hull(){ vector<pt> qs(n * 2); int k = 0; vsort(pts); rep(i, 0, n){ //下側凸包の作成 while(k > 1 && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--; qs[k++] = pts[i]; } for(int i = n - 2, t = k; i >= 0; i--){ //上側凸包の作成 while(k > t && (qs[k-1] - qs[k-2]).det(pts[i] - qs[k-1]) <= 0) k--; qs[k++] = pts[i]; } qs.resize(k-1); n = k-1; pts = qs; return qs; } //点が図形を構成する順番に並んでないと駄目。並んでなければconvex_hullで矯正してから double getArea(){ double s = 0; rep(i, 0, n-1) s += pts[i].det(pts[i+1]); s += pts[n-1].det(pts[0]); return fabs(s) * 0.5; } //http://www004.upp.so-net.ne.jp/s_honma/urawaza/area2.htm double _getArea(){ double s = 0; rep(i, 0, n-1) s += (pts[i].x - pts[i+1].x) * (pts[i].y + pts[i+1].y); s += (pts[n-1].x - pts[0].x) *( pts[n-1].y + pts[0].y); return fabs(s) * 0.5; } //http://imagingsolution.net/math/calc_n_point_area/ //点が時計回りか反時計回りに並んでないと駄目 double getArea2(){ double s = 0; rep(i, 0, n-1) s += pts[i].x * pts[i+1].y + pts[i+1].x * pts[i].y; s += pts[n-1].x * pts[0].y + pts[0].x * pts[n-1].y; return fabs(s) * 0.5; } };
Polygonクラスは座標を表すクラスptのベクタを渡すとそれらの点で結んだ多角形を表します。メンバのptsは座標のベクタ、nは座標数です。
getAreaでその多角形の面積を求めます。この際、座標の順番が辺のつながりの順序に一致していないと正しくない答えが出る場合があります。この場合はptsを凸包にする関数convex_hullを呼ぶと座標が並び替えられます（本来の用途と違うけど）。
PolygonクラスでAOJの問題を解いてみました。
Enclose Pins with a Rubber Band
昔AOJのNo100以下を制覇しようとしてたときに解けなかった問題。凸包を求めればOK。
#include <iostream> #include <cstdio> //#include <cstdlib> #include <cmath> #include <vector> #include <algorithm> #include <string> #include <map> #include <set> #include <queue> //#include <stack> //#include <limits> #include <sstream> //#include <functional> #include <complex> using namespace std; #define len(array) (sizeof (array) / sizeof *(array)) #define rep(i, s, e) for(int i = s;i < e;i++) #define rrep(i, e, s) for(int i = e;s <= i;i--) #define vrange(v) v.begin(), v.end() #define vrrange(v) v.rbegin(), v.rend() #define vsort(v) sort(vrange(v)) #define vrsort(v) sort(vrrange(v)) #define arange(a) a, a + len(a) #define asort(a) sort(arange(a)) #define arsort(a, t) sort(arange(a), greater<t>()) #define afill(a, v) fill(arange(a), v) #define afill2(a, v, t) fill((t *)a, (t *)(a + len(a)), v) #define fmax(a, b) (a < b? b : a) #define fmin(a, b) (a > b? b : a) #define fabs(a) (a < 0? -a : a) #define X real() #define Y imag() typedef unsigned int ui; typedef long long ll; typedef unsigned long long ull; typedef pair<int, int> P; typedef complex<double> p; const int INF = (int)1e9; const int MOD = (int)1e9 + 7; const double EPS = 1e-10; //const int dx[] = {1, -1, 0, 0, 1, -1, -1, 1}; //const int dy[] = {0, 0, 1, -1, -1, -1, 1, 1}; //const int weight[] = {0,1,10,100,1000,10000,100000,1000000,10000000}; //priority_queue< int, vector<int>, greater<int> > q; //二次元座標(from 蟻本) typedef struct _PT{ ... } pt; class Polygon{ ... }; //http://eternuement.blogspot.jp/2011/04/c-string-double-int.html double s2f(string str){ double rt; stringstream ss; ss << str; ss >> rt; return rt; } string f2s(double d){ string rt; stringstream ss; ss << d; ss >> rt; return rt; } vector<string> split( string s, string c ) { vector<string> ret; for(int i = 0, n; i <= s.length(); i = n + 1 ){ n = s.find_first_of( c, i ); if( n == string::npos ) n = s.length(); string tmp = s.substr( i, n-i ); ret.push_back(tmp); } return ret; } void doIt(){ int n; vector<pt> v; string s; vector<string> ss; Polygon pg; while(true){ cin >> n; v.clear(); if(!n) break; rep(i, 0, n){ cin >> s; ss = split(s, ","); v.push_back(pt(s2f(ss[0]), s2f(ss[1]))); } pg = Polygon(v); pg.convex_hull(); cout << n - pg.n << endl; } } int main() { doIt(); return 0; }
Area of Polygon
多角形の面積を求める問題。順序がもとから正しい順序になっているので、座標の順番を変えなくて良い。
#include <iostream> #include <cstdio> //#include <cstdlib> #include <cmath> #include <vector> #include <algorithm> #include <string> #include <map> #include <set> #include <queue> //#include <stack> //#include <limits> #include <sstream> //#include <functional> #include <complex> using namespace std; #define len(array) (sizeof (array) / sizeof *(array)) #define rep(i, s, e) for(int i = s;i < e;i++) #define rrep(i, e, s) for(int i = e;s <= i;i--) #define vrange(v) v.begin(), v.end() #define vrrange(v) v.rbegin(), v.rend() #define vsort(v) sort(vrange(v)) #define vrsort(v) sort(vrrange(v)) #define arange(a) a, a + len(a) #define asort(a) sort(arange(a)) #define arsort(a, t) sort(arange(a), greater<t>()) #define afill(a, v) fill(arange(a), v) #define afill2(a, v, t) fill((t *)a, (t *)(a + len(a)), v) #define fmax(a, b) (a < b? b : a) #define fmin(a, b) (a > b? b : a) #define fabs(a) (a < 0? -a : a) #define X real() #define Y imag() typedef unsigned int ui; typedef long long ll; typedef unsigned long long ull; typedef pair<int, int> P; typedef complex<double> p; const int INF = (int)1e9; const int MOD = (int)1e9 + 7; const double EPS = 1e-10; //const int dx[] = {1, -1, 0, 0, 1, -1, -1, 1}; //const int dy[] = {0, 0, 1, -1, -1, -1, 1, 1}; //const int weight[] = {0,1,10,100,1000,10000,100000,1000000,10000000}; //priority_queue< int, vector<int>, greater<int> > q; //二次元座標(from 蟻本) typedef struct _PT{ ... } pt; class Polygon{ ... }; //http://eternuement.blogspot.jp/2011/04/c-string-double-int.html double s2f(string str){ double rt; stringstream ss; ss << str; ss >> rt; return rt; } string f2s(double d){ string rt; stringstream ss; ss << d; ss >> rt; return rt; } vector<string> split( string s, string c ) { vector<string> ret; for(int i = 0, n; i <= s.length(); i = n + 1 ){ n = s.find_first_of( c, i ); if( n == string::npos ) n = s.length(); string tmp = s.substr( i, n-i ); ret.push_back(tmp); } return ret; } void doIt(){ vector<pt> v; string s; vector<string> ss; Polygon pg; v.clear(); while(!cin.eof()){ cin >> s; ss = split(s, ","); v.push_back(pt(s2f(ss[0]), s2f(ss[1]))); } pg = Polygon(v); pg.convex_hull(); printf("%.7f\n", pg.getArea()); } int main() { doIt(); return 0; }
ミラー・ラビン素数判定法を使ってみた
CodeChefのMay Challenge 2013に挑んだ時に解けなかったWitua and Mathを解いてみました。
Witua and Math
整数Nが与えられる（2 <= N <= 10^18）。オイラーのtotient関数をφとするとき、φ(i)/iを最大とするi(2 <= i <= N)を求めよ。
totient関数のWikipediaを見ると、この問題は「N以下で最大の素数は何か？」という問題と等価である事が分かります。コンテスト中はここまでは分かったのですが、Nが10^18と大きすぎて、制限時間内に可能な素数判定が分かりませんでした。
Editorialによると、ミラー・ロビン素数判定法というアルゴリズムを使えばよいようです。
ただ、このアルゴリズムは確率的素数判定法なので、素数は素数と判定されますが、合成数が素数と判定されてしまうことがあります。
しかし、この配列（Aとおく）のA[i]までの整数は、素数判定法で使うテストケースにi番目の素数までの全ての素数を使う事で決定的に判断できるらしいです（理由は理解できてません）。
この問題では10^18まで判定できればいいので、A[8] < 10^18 < A[9]より、{2,3,5,7,11,13,17,19,23}の9個の素数でテストすればOKです。
Spaghetti Sourceの力を借りて実装。
この問題では乗算によって64ビットをオーバーフローすることがあるので、乗算はEditorialにあるように自前で実装します。
#include <iostream> #include <cstdio> //#include <cstdlib> #include <cmath> #include <vector> #include <algorithm> #include <string> #include <map> #include <set> #include <queue> //#include <stack> //#include <limits> #include <sstream> //#include <functional> #include <complex> using namespace std; #define len(array) (sizeof (array) / sizeof *(array)) #define rep(i, s, e) for(int i = s;i < e;i++) #define rrep(i, e, s) for(int i = e;s <= i;i--) #define vrange(v) v.begin(), v.end() #define vrrange(v) v.rbegin(), v.rend() #define vsort(v) sort(vrange(v)) #define vrsort(v) sort(vrrange(v)) #define arange(a) a, a + len(a) #define asort(a) sort(arange(a)) #define arsort(a, t) sort(arange(a), greater<t>()) #define afill(a, v) fill(arange(a), v) #define afill2(a, v, t) fill((t *)a, (t *)(a + len(a)), v) #define fmax(a, b) (a < b? b : a) #define fmin(a, b) (a > b? b : a) #define fabs(a) (a < 0? -a : a) #define X real() #define Y imag() typedef unsigned int ui; typedef long long ll; typedef unsigned long long ull; typedef pair<int, int> P; typedef complex<double> p; const int INF = (int)1e9; const int MOD = (int)1e9 + 7; const double EPS = 1e-10; //const int dx[] = {1, -1, 0, 0, 1, -1, -1, 1}; //const int dy[] = {0, 0, 1, -1, -1, -1, 1, 1}; //const int weight[] = {0,1,10,100,1000,10000,100000,1000000,10000000}; //priority_queue< int, vector<int>, greater<int> > q; // //calculate (a*b)%m //http://discuss.codechef.com/questions/9723/witmath-editorial ull bigMul(ull a, ull b, ull m){ int base = (int)1e9; ull a_low = a % base, a_high = a / base, b_low = b % base, b_high = b / base, result; result = (a_high * b_high) % m; rep(i, 0, 9) result = (result * 10) % m; result = (result + a_low*b_high + b_low*a_high) % m; rep(i, 0, 9) result = (result * 10) % m; result = (result + a_low*b_low) % m; return result; } //n**p % m ull bigPowMod(ull n, ull p, ull m){ ull ans = 1, ln = n; if(p <= 0) return 1; while(p != 0){ if((p & 1) == 1) ans = bigMul(ans, ln, m); //ans = (ans*ln) % m; //ln = (ln * ln) % m; ln = bigMul(ln, ln, m); p = p >> 1; } return ans; } bool suspect(int a, int s, ull d, ull n) { ull x = bigPowMod(a, d, n); if (x == 1) return true; for (int r = 0; r < s; ++r) { if (x == n - 1) return true; //x = x * x % n; x = bigMul(x, x, n); } return false; } //MillerRabin primality test // {2,7,61,-1} is for n < 4759123141 (= 2^32) // {2,3,5,7,11,13,17,19,23,-1} is for n < 10^16 (at least) int test[] = {2,3,5,7,11,13,17,19,23,-1}; bool MillerRabin(ull n) { if (n <= 1 || (n > 2 && n % 2 == 0)) return false; ull d = n - 1; int s = 0; while (d % 2 == 0) ++s, d /= 2; for (int i = 0; test[i] < n && test[i] != -1; ++i) if (!suspect(test[i], s, d, n)) return false; return true; } void doIt(){ int t; ull n; cin >> t; while(t--){ cin >> n; while(true){ if(MillerRabin(n)){ cout << n << endl; break; } n--; } } } int main() { doIt(); return 0; }
Vimのpegjs用シンタックスファイルを改造してみた
研究でPEG（解析表現文法）を扱っており、パーザージェネレータにPEG.jsを使っています。
pegjs形式のファイルはJavaScriptとPEGの文法規則を並べたものから構成されているので、上手いカラーシンタックスが出来ません。一番近いのはJavaScriptのシンタックスですが、PEGの文法規則の部分は上手く表示してくれません。
そこで、普段使っているEmacs用のpegjs用シンタックスファイルをWebで探してみるも、見つからず。
Vim用とSublime Text2用のシンタックスファイルはあったので、これを機にMacVimを導入してみました。
MacVimはEmacsに比べて全体的にUIが綺麗な印象です。起動もEmacsよりかなり早い。あとEmacsと違ってCtrl-f4を吸収しないのもいいですね。
ただ、やはりモード切り替えが慣れない。ずっと使ってれば慣れるのかな。
こちらのpegjs用シンタックスを使ってみましたが、JavaScriptの部分のコードの色づけをしてくれません。僕が普段書くコードはJavaScriptの部分が大半を占めるので、これではあまり使い勝手がよくありません。
そこで、マニュアルを見ながらpegjs.vimファイルを改造してJavaScriptの部分にも色が付くようにしました。
具体的には、MacVimに標準で入っているJavaScript用シンタックスファイルのjavascript.vimを丸ごとpegjs.vimにコピーし、pegjs.vimでJavaScriptと判断される部分に適用されるようにしました。
"This file is pegjs.vim + javascript.vim(vim default). " pegjs.vim " Language: PEG JS Grammars " Maintainer: Andrew Lunny <alunny@gmail.com> " Latest Revision: 15 June 2012 " License: MIT " PEG's rule Syntax(import from pegjs.vim) if exists("b:current_syntax") finish endif syn match pegOperator "[/.+*?&!]" syn region charSet start="\[" end="\]" syn match rule "[a-zA-Z$_][a-zA-Z$_0-9]*" syn match ruleDef "[_a-zA-Z$][a-zA-Z$_0-9]*" contained syn match equals "=" contained syn match initialize "[a-zA-Z$_][a-zA-Z$_0-9]*[\n\t ]*\".*\"[\n\t ]*=" contains=ruleDef,equals,innerLiteral syn match initialize "[a-zA-Z$_][a-zA-Z$_0-9]*[\n\t ]*'.*'[\n\t ]*=" contains=ruleDef,equals,innerLiteral syn match initialize "[a-zA-Z$_][a-zA-Z$_0-9]*[\n\t ]*=" contains=ruleDef,equals syn match exprLabel "[a-zA-Z$_][a-zA-Z$_0-9]*:"he=e-1 syn region literal start="'" end="'" syn region literal start="\"" end="\"" syn region innerLiteral start="'" end="'" contained syn region innerLiteral start="\"" end="\"" contained syn region comment start="/[*]" end="[*]/" syn region comment start="//" end="\n" syn region jsBlock start="{" end="}" contains=ALLBUT,pegOperator,charSet,rule,ruleDef,equals,initialize,exprLabel,literal,innerLiteral,comment,jsBlock hi def link ruleDef PreProc hi def link rule Type hi def link exprLabel Identifier hi def link literal String hi def link pegOperator Operator hi def link charSet Operator hi def link whitespace PreProc hi def link innerLiteral Comment "hi def link jsBlock Comment "hi jsBlock guibg=#cccccc "hi ruleDef gui=underline,bold " JavaScript Block Syntax(import from javascript.vim) if !exists("main_syntax") if version < 600 syntax clear elseif exists("b:current_syntax") finish endif let main_syntax = 'javascript' endif " Drop fold if it set but vim doesn't support it. if version < 600 && exists("javaScript_fold") unlet javaScript_fold endif syn keyword javaScriptCommentTodo TODO FIXME XXX TBD contained syn match javaScriptLineComment "\/\/.*" contains=@Spell,javaScriptCommentTodo contained syn match javaScriptCommentSkip "^[ \t]*\*\($\|[ \t]\+\)" contained syn region javaScriptComment start="/\*" end="\*/" contains=@Spell,javaScriptCommentTodo contained syn match javaScriptSpecial "\\\d\d\d\|\\." contained syn region javaScriptStringD start=+"+ skip=+\\\\\|\\"+ end=+"\|$+ contains=javaScriptSpecial,@htmlPreproc contained syn region javaScriptStringS start=+'+ skip=+\\\\\|\\'+ end=+'\|$+ contains=javaScriptSpecial,@htmlPreproc contained syn match javaScriptSpecialCharacter "'\\.'" contained syn match javaScriptNumber "-\=\<\d\+L\=\>\|0[xX][0-9a-fA-F]\+\>" syn region javaScriptRegexpString start=+/[^/*]+me=e-1 skip=+\\\\\|\\/+ end=+/[gi]\{0,2\}\s*$+ end=+/[gi]\{0,2\}\s*[;.,)\]}]+me=e-1 contains=@htmlPreproc oneline contained syn keyword javaScriptConditional if else switch syn keyword javaScriptRepeat while for do in syn keyword javaScriptBranch break continue syn keyword javaScriptOperator new delete instanceof typeof syn keyword javaScriptType Array Boolean Date Function Number Object String RegExp syn keyword javaScriptStatement return with syn keyword javaScriptBoolean true false syn keyword javaScriptNull null undefined syn keyword javaScriptIdentifier arguments this var let syn keyword javaScriptLabel case default syn keyword javaScriptException try catch finally throw syn keyword javaScriptMessage alert confirm prompt status syn keyword javaScriptGlobal self window top parent syn keyword javaScriptMember document event location syn keyword javaScriptDeprecated escape unescape syn keyword javaScriptReserved abstract boolean byte char class const debugger double enum export extends final float goto implements import int interface long native package private protected public short static super synchronized throws transient volatile if exists("javaScript_fold") syn match javaScriptFunction "\<function\>" syn region javaScriptFunctionFold start="\<function\>.*[^};]$" end="^\z1}.*$" transparent fold keepend syn sync match javaScriptSync grouphere javaScriptFunctionFold "\<function\>" syn sync match javaScriptSync grouphere NONE "^}" setlocal foldmethod=syntax setlocal foldtext=getline(v:foldstart) else syn keyword javaScriptFunction function "syn match javaScriptBraces "[{}\[\]]" contained syn match javaScriptBraces "[\[\]]" contained syn region jsBlocks matchgroup=Function start="{" end="}" contains=ALLBUT,pegOperator,charSet,rule,ruleDef,equals,initialize,exprLabel,literal,innerLiteral,comment,jsBlock containedin=jsBlock,jsBlocks syn match javaScriptParens "[()]" endif syn sync fromstart syn sync maxlines=100 if main_syntax == "javascript" syn sync ccomment javaScriptComment endif " Define the default highlighting. " For version 5.7 and earlier: only when not done already " For version 5.8 and later: only when an item doesn't have highlighting yet if version >= 508 || !exists("did_javascript_syn_inits") if version < 508 let did_javascript_syn_inits = 1 command -nargs=+ HiLink hi link <args> else command -nargs=+ HiLink hi def link <args> endif HiLink javaScriptComment Comment HiLink javaScriptLineComment Comment HiLink javaScriptCommentTodo Todo HiLink javaScriptSpecial Special HiLink javaScriptStringS String HiLink javaScriptStringD String HiLink javaScriptCharacter Character HiLink javaScriptSpecialCharacter javaScriptSpecial HiLink javaScriptNumber javaScriptValue HiLink javaScriptConditional Conditional HiLink javaScriptRepeat Repeat HiLink javaScriptBranch Conditional HiLink javaScriptOperator Operator HiLink javaScriptType Type HiLink javaScriptStatement Statement HiLink javaScriptFunction Function HiLink javaScriptBraces Function HiLink javaScriptError Error HiLink javaScrParenError javaScriptError HiLink javaScriptNull Keyword HiLink javaScriptBoolean Boolean HiLink javaScriptRegexpString String HiLink javaScriptIdentifier Identifier HiLink javaScriptLabel Label HiLink javaScriptException Exception HiLink javaScriptMessage Keyword HiLink javaScriptGlobal Keyword HiLink javaScriptMember Keyword HiLink javaScriptDeprecated Exception HiLink javaScriptReserved Keyword HiLink javaScriptDebug Debug HiLink javaScriptConstant Label delcommand HiLink endif let b:current_syntax = "javascript" if main_syntax == 'javascript' unlet main_syntax endif
上記を.vim/syntax/pegjs.vimに書いて、.vimrcに
autocmd BufNewFile,BufReadPost *.pegjs set filetype=pegjs
と書くと拡張子がpegjsのファイルにこのシンタックスが適用されます。今のところ問題なく動くみたい。しばらくvimを使ってみます。
偶然にも今日はpegjs.vimの作者様のLatest Rivisionのぴったり1年後でした。
