SRM272 DIV2 1000 vectorpolygon

問題

　N個の2次元ベクトルが与えられる。このベクトルからいくつかを選び適切な順に並べ直し、凸多角形を作る。作ることのできる凸多角形のうち面積が最大のものの面積を求めよ。凸多角形が作れない場合は0を出力せよ。
　1≦N≦8

方針

　x軸となす角度が小さい順にベクトルを並べ替える。適当に選んだベクトル集合の和がゼロベクトルなら閉じた領域を作れる。 　あとは凸判定だけど、x軸となす角が小さい順に辺が並んでいるので凸になっている。（逆に凸多角形ならば、辺を反時計まわりにみれば辺とx軸がなす角はかならず増えていく。） 　面積は外積を使って計算できる。

ソースコード

view plainprint?

1. using namespace std;  
2.   
3. #define REP(i, n) for(int i=0; i<(int)(n); i++)  
4. #define FOR(i, s, e) for (int i = (int)(s); i < (int)(e); i++)  
5.   
6. class VectorPolygon {  
7.     template <class T> T crossProduct(T x1, T y1, T x2, T y2) {  
8.         return x1*y2 - x2*y1;  
9.     }  
10.   
11.     double getArea(vector<int> &x, vector<int> &y) {  
12.         int N = x.size();  
13.   
14.         long long ret = 0;  
15.         REP(i, N) {  
16.             int j = (i+1)%N;  
17.             ret += crossProduct<long long>(x[i], y[i], x[j], y[j]);  
18.         }  
19.   
20.         return 0.5*ret;  
21.     }  
22.   
23.     double getAngle(int x, int y) {  
24.         double ret = atan2(y, x);  
25.         if (ret >= 0)  
26.             return ret;  
27.         else  
28.             return 2*M_PI+ret;  
29.     }  
30.   
31.   
32. public:  
33.     double maxArea(vector<int> dx, vector<int> dy) {  
34.         int N = dx.size();  
35.         REP(i, N) FOR(j, i+1, N) {  
36.             if (getAngle(dx[i], dy[i]) > getAngle(dx[j], dy[j])) {  
37.                 swap(dx[i], dx[j]);  
38.                 swap(dy[i], dy[j]);  
39.             }  
40.         }  
41.   
42.         REP(i, N)  
43.             cout << dx[i] << " " << dy[i] << " " << getAngle(dx[i], dy[i]) << endl;  
44.   
45.         double ret = 0.0;  
46.         vector<int> vx;  
47.         vector<int> vy;  
48.   
49.         REP(mask, 1<<N) {  
50.             vx.clear();  
51.             vy.clear();  
52.             int x = 0;  
53.             int y = 0;  
54.             int z = 0;  
55.   
56.             vx.push_back(0);  
57.             vy.push_back(0);  
58.             REP(i, N) if (mask & 1<<i) {  
59.                 x += dx[i];  
60.                 y += dy[i];  
61.                 if (!x && !y)  
62.                     ++z;  
63.                 else {  
64.                     vx.push_back(x);  
65.                     vy.push_back(y);  
66.                 }  
67.             }  
68.   
69.             if (x || y || (int)vx.size() < 3) continue;  
70.             if (z != 1) continue;  
71.   
72.             ret = max(ret, getArea(vx, vy));  
73.         }  
74.   
75.         return ret;  
76.     }  
77. };  

using namespace std;

#define REP(i, n) for(int i=0; i<(int)(n); i++)
#define FOR(i, s, e) for (int i = (int)(s); i < (int)(e); i++)

class VectorPolygon {
    template <class T> T crossProduct(T x1, T y1, T x2, T y2) {
        return x1*y2 - x2*y1;
    }

    double getArea(vector<int> &x, vector<int> &y) {
        int N = x.size();

        long long ret = 0;
        REP(i, N) {
            int j = (i+1)%N;
            ret += crossProduct<long long>(x[i], y[i], x[j], y[j]);
        }

        return 0.5*ret;
    }

    double getAngle(int x, int y) {
        double ret = atan2(y, x);
        if (ret >= 0)
            return ret;
        else
            return 2*M_PI+ret;
    }


public:
    double maxArea(vector<int> dx, vector<int> dy) {
        int N = dx.size();
        REP(i, N) FOR(j, i+1, N) {
            if (getAngle(dx[i], dy[i]) > getAngle(dx[j], dy[j])) {
                swap(dx[i], dx[j]);
                swap(dy[i], dy[j]);
            }
        }

        REP(i, N)
            cout << dx[i] << " " << dy[i] << " " << getAngle(dx[i], dy[i]) << endl;

        double ret = 0.0;
        vector<int> vx;
        vector<int> vy;

        REP(mask, 1<<N) {
            vx.clear();
            vy.clear();
            int x = 0;
            int y = 0;
            int z = 0;

            vx.push_back(0);
            vy.push_back(0);
            REP(i, N) if (mask & 1<<i) {
                x += dx[i];
                y += dy[i];
                if (!x && !y)
                    ++z;
                else {
                    vx.push_back(x);
                    vy.push_back(y);
                }
            }

            if (x || y || (int)vx.size() < 3) continue;
            if (z != 1) continue;

            ret = max(ret, getArea(vx, vy));
        }

        return ret;
    }
};

読書「1Q84 BOOK1 〈4月‐6月〉前編」

村上春樹の「1Q84 BOOK1 〈4月‐6月〉前編」 を読んだ。村上春樹は友達が好きで何度も勧められていた。英会話のフィリピン人の講師や、アメリカ人のPen Palも好きと言っていたので、これは読むしかないと思ってよんでみた。

ストーリーは良く練られていておもしろく、かつ文章の表現も印象に残るものが多くいろんな人が好きだという理由が分かりました。女殺し屋の青豆と、塾で働く数学講師かつ趣味で小説を書いている天吾の二人に焦点をあてた群像劇のような感じで章立てされています。加えて、美少女不思議系女子高生フカエリや、青豆に殺しを依頼する老婦人などなど個性豊かなキャラクターたちが物語を一層おもしろいものにします。

村上春樹の小説は読みにくいとか理解できないという人が多いですが（そういう人達の大半は「海辺のカフカ」を読んだといいます。）、この「1Q84」はとても読みやすいです。文庫本が出てるので興味のある人は是非。

Facebook Hacker Cup 2011 Round2 Scott' New Trick

問題概要

　素数pが与えられる。サイズnとサイズmの整数列A, Bを考える。但しA, Bは以下のように決まる。
A[0] = A1
A[1] = A2
A[i] = (A3 * A[i-2] + A4 * A[i-1] * A5) % p, i = 2,3,..,n-1
B[0] = B1
B[1] = B2
B[i] = (B3 * A[i-2] + B4 * A[i-1] * B5) % p, i = 2,3,..,m-1

(A[i] * B[j]) mod p < lとなるような(i, j)の組の個数を求めよ。

2 ≤ P < 250,000
1 ≤ L ≤ P
2 ≤ N, M ≤ 10,000,000
0 ≤ A1, A2, A3, A4, A5, B1, B2, B3, B4, B5 < P

方針

　ここのサイトがとても分かりやすい。一応日本語で概要だけまとめておく。

　素数pには必ず原始根が存在する。以下では、素数pの原始根をgと記す。 x * y (mod p) = z となるようなx, yを考えよ。という問題は、x, y, zをgのべき乗で表してやると（つまりg^x' * g^y' = g^z'と表してやると）、 x' + y' (mod p) = z' となるような(x', y')を考えよという問題に読み変えることができる。
　よって、まず原始根を求めて[1, p-1]までの値 -> それをgのべき乗で表した場合べきの数はいくらか？というlookupテーブルを作成する。原始根の計算は愚直にやると、p-2個の候補のp-1個のべき乗を調べることになるので、O(p²)くらいかかる（ただし、2が原始根になる場合もあるので、Ω(p)）。エラトステネスの篩 + 位数がφ(p) = p-1の約数になること + 繰り返し二乗法を利用すれば、(p-1)を素因数で割ったものだけを調べることにより、O(p * log²(p))くらいで計算できそうだったが、予め最小の原始根を調べてみると、せいぜい40程度だったので、愚直な実装のままでいくことにした。
　次に、実際にA, Bを生成しながら、どのべき数がどれだけ数列内に含まれているかという情報を保存していく。Aに含まれるべき数がiのものの個数をw[i], Bに含まれるべき数がjのものの個数をv[j]とする。 このとき、Aの要素とBの要素を掛け合わせた結果のべき数がiとなるようなもの組み合わせの数は、

のように表すことができる。これは畳みこみ和の形をしているので、カラツバ法もしくはFFTを用いて高速に計算することができる。今回はFFTで計算した。あとは、gⁱがlより小さいかどうか判定して、小さければu[i]を答えに足していくだけ。また、ゼロは原始根のべき乗で表せないので例外的に処理をした。

原始根って何？原始根は知ってるけど、何ですべての素数に原始根が存在するの？という人には、以下の本がお勧めです。

ソースコード

　ジャッジ用のテストデータがダウンロードできなかったため、完全な正誤判定はできていない。（ランダムに比較的小さなサイズのテストケースを生成して愚直な解法と比べて検算した。また、大きなテストケースも生成して許容時間内（6分内）にすべてのテストケースを計算できることを確認した。）誰かローカルにテストデータ持ってたら、ください。。

view plainprint?

1. #include <iostream>  
2. #include <vector>  
3. #include <algorithm>  
4. #include <complex>  
5.   
6. using namespace std;  
7.   
8. /* 
9.  * get a primitive root 
10.  */  
11. int getPrimitiveRoot(int p) {  
12.     for (int i = 1; i < p; i++) {  
13.         long long x = i;  
14.         bool ck = true;  
15.         for (int j = 1; j < p - 1; j++) {  
16.             if (x == 1) {  
17.                 ck = false;  
18.                 break;  
19.             }  
20.             x = x * i % p;  
21.         }  
22.         if (ck)  
23.             return i;  
24.     }  
25.     return -1;  
26. }  
27.   
28. /* 
29.  * Fast Fourier Transform 
30.  */  
31. typedef complex<long double> Complex;  
32. const int TIME_TO_FREQ = -1;  
33. const int FREQ_TO_TIME = 1;  
34.   
35. void fft(Complex x[], int N, int sgn) {  
36.     Complex theta = Complex(0, sgn * 2 * M_PI / N);  
37.   
38.     for (int m = N; m >= 2; m >>= 1) {  
39.         int mh = m >> 1;  
40.         for (int i = 0; i < mh; i++) {  
41.             Complex w = exp((long double)i * theta);  
42.             for (int j = i; j < N; j += m) {  
43.                 int k = j + mh;  
44.                 Complex tmp = x[j] - x[k];  
45.                 x[j] += x[k];  
46.                 x[k] = w * tmp;  
47.             }  
48.         }  
49.         theta *= 2;  
50.     }  
51.   
52.     int i = 0;  
53.     for (int j = 1; j < N - 1; j++) {  
54.         for (int k = N >> 1; k > (i ^= k); k >>= 1) ;  
55.         if (j < i) swap(x[i], x[j]);  
56.     }  
57.   
58.     if (sgn == FREQ_TO_TIME) for (int i = 0; i < N; i++) x[i] /= N;  
59. }  
60.   
61. /* 
62.  * main solver 
63.  */  
64. const int MAX_SEQ_SIZE = 10000000;  
65. const int MAX_PRIME_NUM = 1<<18;  
66. long long a[MAX_SEQ_SIZE];  
67. long long b[MAX_SEQ_SIZE];  
68. int discreteLog[MAX_PRIME_NUM];  
69. int powVal[MAX_PRIME_NUM];  
70. Complex w[2*MAX_PRIME_NUM];  
71. Complex v[2*MAX_PRIME_NUM];  
72. Complex u[2*MAX_PRIME_NUM];  
73.   
74. long long solve(int p, int l, int n, int m) {  
75.     // calculate a primitive root on modular p  
76.     // and make a look-up table (n, discrete log(n))  
77.     int g = getPrimitiveRoot(p);  
78.     long long x = 1;  
79.     for (int i = 0; i < p - 1; i++) {  
80.         discreteLog[x] = i;  
81.         powVal[i] = x;  
82.         x = x * g % p;  
83.     }  
84.   
85.     // calculate convolutions with FFTs  
86.     int N = 1;  
87.     while (N <= p)  
88.         N *= 2;  
89.     N *= 2;  
90.   
91.     for (int i = 0; i < N; i++) {  
92.         w[i] = 0;  
93.         v[i] = 0;  
94.     }  
95.   
96.     long long ret = 0;  
97.     long long z1 = 0;  
98.     long long z2 = 0;  
99.     for (int i = 0; i < n; i++) {  
100.         if (!a[i]) {  
101.             ret += m;  
102.             ++z1;  
103.         }  
104.         else  
105.             w[discreteLog[a[i]]].real()++;  
106.     }  
107.     for (int i = 0; i < m; i++) {  
108.         if (!b[i]) {  
109.             ret += n;  
110.             ++z2;  
111.         }  
112.         else  
113.             v[discreteLog[b[i]]].real()++;  
114.     }  
115.     ret -= z1 * z2;  
116.   
117.     fft(w, N, TIME_TO_FREQ);  
118.     fft(v, N, TIME_TO_FREQ);  
119.   
120.     for (int i = 0; i < N; i++)  
121.         u[i] = w[i] * v[i];  
122.   
123.     fft(u, N, FREQ_TO_TIME);  
124.   
125.     for (int i = 0; i < p-1; i++) {  
126.         if (powVal[i] >= l)  
127.             continue;  
128.         if (i == p-2)  
129.             ret += (long long)(u[i].real() + 0.5);  
130.         else  
131.             ret += (long long)(u[i].real() + u[i+p-1].real() + 0.5);  
132.     }  
133.   
134.     return ret;  
135. }  
136.   
137. /* 
138.  * main (for input and output only) 
139.  */  
140. #include <float.h>  
141. int main() {  
142.     int T;  
143.     cin >> T;  
144.   
145.     for (int t = 0; t < T; t++) {  
146.         cerr << "Case #" << t+1 << " solving....." << endl;  
147.         int p, l;  
148.         cin >> p >> l;  
149.   
150.         // generate sequence A  
151.         int n, a3, a4, a5;  
152.         cin >> n >> a[0] >> a[1] >> a3 >> a4 >> a5;  
153.         for (int i = 2; i < n; i++)  
154.             a[i] = (a[i - 2] * a3 + a[i - 1] * a4 + a5) % p;  
155.   
156.         // generate sequence B  
157.         int m, b3, b4, b5;  
158.         cin >> m >> b[0] >> b[1] >> b3 >> b4 >> b5;  
159.         for (int i = 2; i < m; i++)  
160.             b[i] = (b[i - 2] * b3 + b[i - 1] * b4 + b5) % p;  
161.   
162.         long long ret = solve(p, l, n, m);  
163.         cout << "Case #" << t+1 << ": " << ret << endl;  
164.     }  
165.   
166.     return 0;  
167. }  

#include <iostream>
#include <vector>
#include <algorithm>
#include <complex>

using namespace std;

/*
 * get a primitive root
 */
int getPrimitiveRoot(int p) {
    for (int i = 1; i < p; i++) {
        long long x = i;
        bool ck = true;
        for (int j = 1; j < p - 1; j++) {
            if (x == 1) {
                ck = false;
                break;
            }
            x = x * i % p;
        }
        if (ck)
            return i;
    }
    return -1;
}

/*
 * Fast Fourier Transform
 */
typedef complex<long double> Complex;
const int TIME_TO_FREQ = -1;
const int FREQ_TO_TIME = 1;

void fft(Complex x[], int N, int sgn) {
    Complex theta = Complex(0, sgn * 2 * M_PI / N);

    for (int m = N; m >= 2; m >>= 1) {
        int mh = m >> 1;
        for (int i = 0; i < mh; i++) {
            Complex w = exp((long double)i * theta);
            for (int j = i; j < N; j += m) {
                int k = j + mh;
                Complex tmp = x[j] - x[k];
                x[j] += x[k];
                x[k] = w * tmp;
            }
        }
        theta *= 2;
    }

    int i = 0;
    for (int j = 1; j < N - 1; j++) {
        for (int k = N >> 1; k > (i ^= k); k >>= 1) ;
        if (j < i) swap(x[i], x[j]);
    }

    if (sgn == FREQ_TO_TIME) for (int i = 0; i < N; i++) x[i] /= N;
}

/*
 * main solver
 */
const int MAX_SEQ_SIZE = 10000000;
const int MAX_PRIME_NUM = 1<<18;
long long a[MAX_SEQ_SIZE];
long long b[MAX_SEQ_SIZE];
int discreteLog[MAX_PRIME_NUM];
int powVal[MAX_PRIME_NUM];
Complex w[2*MAX_PRIME_NUM];
Complex v[2*MAX_PRIME_NUM];
Complex u[2*MAX_PRIME_NUM];

long long solve(int p, int l, int n, int m) {
    // calculate a primitive root on modular p
    // and make a look-up table (n, discrete log(n))
    int g = getPrimitiveRoot(p);
    long long x = 1;
    for (int i = 0; i < p - 1; i++) {
        discreteLog[x] = i;
        powVal[i] = x;
        x = x * g % p;
    }

    // calculate convolutions with FFTs
    int N = 1;
    while (N <= p)
        N *= 2;
    N *= 2;

    for (int i = 0; i < N; i++) {
        w[i] = 0;
        v[i] = 0;
    }

    long long ret = 0;
    long long z1 = 0;
    long long z2 = 0;
    for (int i = 0; i < n; i++) {
        if (!a[i]) {
            ret += m;
            ++z1;
        }
        else
            w[discreteLog[a[i]]].real()++;
    }
    for (int i = 0; i < m; i++) {
        if (!b[i]) {
            ret += n;
            ++z2;
        }
        else
            v[discreteLog[b[i]]].real()++;
    }
    ret -= z1 * z2;

    fft(w, N, TIME_TO_FREQ);
    fft(v, N, TIME_TO_FREQ);

    for (int i = 0; i < N; i++)
        u[i] = w[i] * v[i];

    fft(u, N, FREQ_TO_TIME);

    for (int i = 0; i < p-1; i++) {
        if (powVal[i] >= l)
            continue;
        if (i == p-2)
            ret += (long long)(u[i].real() + 0.5);
        else
            ret += (long long)(u[i].real() + u[i+p-1].real() + 0.5);
    }

    return ret;
}

/*
 * main (for input and output only)
 */
#include <float.h>
int main() {
    int T;
    cin >> T;

    for (int t = 0; t < T; t++) {
        cerr << "Case #" << t+1 << " solving....." << endl;
        int p, l;
        cin >> p >> l;

        // generate sequence A
        int n, a3, a4, a5;
        cin >> n >> a[0] >> a[1] >> a3 >> a4 >> a5;
        for (int i = 2; i < n; i++)
            a[i] = (a[i - 2] * a3 + a[i - 1] * a4 + a5) % p;

        // generate sequence B
        int m, b3, b4, b5;
        cin >> m >> b[0] >> b[1] >> b3 >> b4 >> b5;
        for (int i = 2; i < m; i++)
            b[i] = (b[i - 2] * b3 + b[i - 1] * b4 + b5) % p;

        long long ret = solve(p, l, n, m);
        cout << "Case #" << t+1 << ": " << ret << endl;
    }

    return 0;
}

ICPC2006 アジア地区予選 Enjoyable Commutation

問題概要

　グラフ(V, E)が与えられる。a ∈ Vからb∈Vまでの経路のうちk番目に短い経路を生成せよ。ただし、1つの経路内で同じ節点を2回以上通ってはいけない。同じ長さの経路が2つ以上ある場合は、通過する節点を順番に並べたときに辞書順最小のものが一番短い経路であると判断する。

　2 ≦ n ≦ 50
　0 ≦ m ≦n(n-1)
　1 ≦ k ≦ 200

方針

　k-th shortest pathという典型問題。基本的にはダイクストラの要領でパスを生成していって、A*探索みたいな感じで最終的にえられるパスの長さの下限が小さいものから吸い上げていけばいい。
　ややこしいのが、同じ節点を複数回通ってはいけないという制約。ビットマスクで訪れた点を持っておきパス生成のときまだ訪れていない節点のみに遷移するのは当然やるとして、これに加えてそこからゴールにたどり着くまでの下限値も変化することに注意。よって、毎回priority_queueからpopされた状態におけるグラフ形状でダイクストラをしてゴールまでのパスの長さの下限値を出すというようなことをやってみた。グラフの形状が変化しなければ（同じ節点を何度通ってもいいならば）、最初にゴールから逆方向にダイクストラをすれば経路の下限は出るが、今回の場合はそのやり方だと精度の悪い下限値が出てしまい現実的な時間で計算が完了しないケースがあった。A*は評価関数の精度がいいほど速く計算できるが、今回もそれと同じだなと思った。
　あと、プラスα的な高速化として、すでにk回popされた節点にはもうpushしないようにするということをした。

※追記2012/04/19 この問題がPOJ(problem 3137)に登録されているのを発見しました。submitしたらTLEでした。さらなる高速化が必要なようです。

ソースコード

　スパゲッティなソースを載せておきます。

view plainprint?

1. #include <iostream>  
2. #include <vector>  
3. #include <queue>  
4. #include <sstream>  
5. #include <cstring>  
6.   
7. using namespace std;  
8.   
9. struct edge {  
10.     int to;  
11.     int cost;  
12. };  
13.   
14. typedef pair<int, int> PII;  
15. const int MAX_V = 50;  
16. const int INF = (int)1e9;  
17.   
18. int V;  
19. vector<edge> G[MAX_V];  
20. int dist[MAX_V];  
21.   
22. /* 
23.  * s: starting point 
24.  * t: ending point 
25.  */  
26. void dijkstra(int s, int t, long long vis) {  
27.     priority_queue<PII, vector<PII>, greater<PII> > Q;  
28.     fill(dist, dist+V, INF);  
29.     dist[s] = 0;  
30.     Q.push(PII(0, s));  
31.   
32.     while(!Q.empty()) {  
33.         int d = Q.top().first;  
34.         int v = Q.top().second;  
35.         Q.pop();  
36.   
37.         if (d > dist[v]) continue;  
38.         if (v == t) return;  
39.   
40.         for (int i = 0; i < (int)G[v].size(); i++) {  
41.             int to = G[v][i].to;  
42.             int cost = G[v][i].cost;  
43.   
44.             if (d + cost < dist[to] && !(vis & 1LL << to)) {  
45.                 dist[to] = d + cost;  
46.                 Q.push(PII(dist[to], to));  
47.             }  
48.         }  
49.     }  
50. }  
51.   
52. class Status {  
53. public:  
54.     int v;  
55.     int len;  
56.     int lb;  
57.     long long vis;  
58.     vector<int> path;  
59.   
60.     Status(int v, int len, long long vis, vector<int>path, int goal) {  
61.         this->v = v;  
62.         this->len = len;  
63.         this->vis = vis;  
64.         dijkstra(v, goal, vis);  
65.         this->lb = min(INF, len + dist[goal]);  
66.         this->path = path;  
67.     }  
68.   
69.     bool operator>(const Status &v) const {  
70.         if (this->lb != v.lb)  
71.             return this->lb > v.lb;  
72.         return this->path > v.path;  
73.     }  
74. };  
75.   
76. int cntVis[MAX_V];  
77. vector<int> kthShortestPath(int s, int t, int k) {  
78.     priority_queue<Status, vector<Status>, greater<Status> > Q;  
79.     vector<int> path;  
80.     path.push_back(s);  
81.     Q.push(Status(s, 0, 1LL<<s, path, t));  
82.     memset(cntVis, 0, sizeof(cntVis));  
83.     int m = 0;  
84.   
85.     while (!Q.empty()) {  
86.         int v = Q.top().v;  
87.         int len = Q.top().len;  
88.         long long vis = Q.top().vis;  
89.         vector<int>path = Q.top().path;  
90.         Q.pop();  
91.   
92.         ++cntVis[v];  
93.         if (v == t) {  
94.             if (++m == k)  
95.                 return path;  
96.             continue;  
97.         }  
98.   
99.         path.push_back(-1);  
100.         for (int i = 0; i < (int)G[v].size(); i++) {  
101.             int to = G[v][i].to;  
102.             int cost = G[v][i].cost;  
103.   
104.             if (!(vis & 1LL<<to) && cntVis[to] < k) {  
105.                 path.back() = to;  
106.                 Q.push(Status(to, len+cost, vis|1LL<<to, path, t));  
107.             }  
108.         }  
109.     }  
110.   
111.     return vector<int>();  
112. }  
113.   
114. int main() {  
115.     int n, m, k, a, b;  
116.     int x, y, d;  
117.   
118.     while (cin >> n >> m >> k >> a >> b) {  
119.         if (n+m+k+a+b == 0)  
120.             break;  
121.   
122.         --a;  
123.         --b;  
124.         V = n;  
125.   
126.         // construct a graph  
127.         for (int i = 0; i < V; i++)  
128.             G[i].clear();  
129.         for (int i = 0; i < m; i++) {  
130.             cin >> x >> y >> d;  
131.             G[--x].push_back((edge){--y, d});  
132.         }  
133.   
134.         vector<int> ret = kthShortestPath(a, b, k);  
135.         if (ret.empty())  
136.             cout << "None" << endl;  
137.         else {  
138.             stringstream ss;  
139.             for (int i = 0; i < (int)ret.size(); i++) {  
140.                 if (i) ss << "-";  
141.                 ss << ret[i]+1;  
142.             }  
143.             cout << ss.str() << endl;  
144.         }  
145.     }  
146.   
147.     return 0;  
148. }  

#include <iostream>
#include <vector>
#include <queue>
#include <sstream>
#include <cstring>

using namespace std;

struct edge {
    int to;
    int cost;
};

typedef pair<int, int> PII;
const int MAX_V = 50;
const int INF = (int)1e9;

int V;
vector<edge> G[MAX_V];
int dist[MAX_V];

/*
 * s: starting point
 * t: ending point
 */
void dijkstra(int s, int t, long long vis) {
    priority_queue<PII, vector<PII>, greater<PII> > Q;
    fill(dist, dist+V, INF);
    dist[s] = 0;
    Q.push(PII(0, s));

    while(!Q.empty()) {
        int d = Q.top().first;
        int v = Q.top().second;
        Q.pop();

        if (d > dist[v]) continue;
        if (v == t) return;

        for (int i = 0; i < (int)G[v].size(); i++) {
            int to = G[v][i].to;
            int cost = G[v][i].cost;

            if (d + cost < dist[to] && !(vis & 1LL << to)) {
                dist[to] = d + cost;
                Q.push(PII(dist[to], to));
            }
        }
    }
}

class Status {
public:
    int v;
    int len;
    int lb;
    long long vis;
    vector<int> path;

    Status(int v, int len, long long vis, vector<int>path, int goal) {
        this->v = v;
        this->len = len;
        this->vis = vis;
        dijkstra(v, goal, vis);
        this->lb = min(INF, len + dist[goal]);
        this->path = path;
    }

    bool operator>(const Status &v) const {
        if (this->lb != v.lb)
            return this->lb > v.lb;
        return this->path > v.path;
    }
};

int cntVis[MAX_V];
vector<int> kthShortestPath(int s, int t, int k) {
    priority_queue<Status, vector<Status>, greater<Status> > Q;
    vector<int> path;
    path.push_back(s);
    Q.push(Status(s, 0, 1LL<<s, path, t));
    memset(cntVis, 0, sizeof(cntVis));
    int m = 0;

    while (!Q.empty()) {
        int v = Q.top().v;
        int len = Q.top().len;
        long long vis = Q.top().vis;
        vector<int>path = Q.top().path;
        Q.pop();

        ++cntVis[v];
        if (v == t) {
            if (++m == k)
                return path;
            continue;
        }

        path.push_back(-1);
        for (int i = 0; i < (int)G[v].size(); i++) {
            int to = G[v][i].to;
            int cost = G[v][i].cost;

            if (!(vis & 1LL<<to) && cntVis[to] < k) {
                path.back() = to;
                Q.push(Status(to, len+cost, vis|1LL<<to, path, t));
            }
        }
    }

    return vector<int>();
}

int main() {
    int n, m, k, a, b;
    int x, y, d;

    while (cin >> n >> m >> k >> a >> b) {
        if (n+m+k+a+b == 0)
            break;

        --a;
        --b;
        V = n;

        // construct a graph
        for (int i = 0; i < V; i++)
            G[i].clear();
        for (int i = 0; i < m; i++) {
            cin >> x >> y >> d;
            G[--x].push_back((edge){--y, d});
        }

        vector<int> ret = kthShortestPath(a, b, k);
        if (ret.empty())
            cout << "None" << endl;
        else {
            stringstream ss;
            for (int i = 0; i < (int)ret.size(); i++) {
                if (i) ss << "-";
                ss << ret[i]+1;
            }
            cout << ss.str() << endl;
        }
    }

    return 0;
}

読書「ゴールデンスランバー」

伊坂幸太郎のゴールデンスランバーを読みました。J・F・ケネディ暗殺事件からヒントを得たようなフィクション小説です。日本国の首相が爆弾で暗殺されて何故か自分が容疑者として指名手配されます。次々と出てくる証言者、行ったことのない店の防犯カメラに何故か映っている自分の姿、異常なほどに強引な警察の捜査。
　
　最後の締め方が面白かったです。「痴漢は死ね。」「ご主人浮気してますよ。」「たいへんよくできました。」なんというか味があって思わずにんまりしました。ちなみにタイトルはビートルズのアビーロードの中の曲名から取ったそうです。はて？そんな曲あったかな。と思って改めて聞いてみると名曲でした。

SRM 283 DIV2 600 PowerSupply

問題

　デカルト平面上にN個の点(x_i, y_i)がある。そこにx-軸、y-軸、2つの対角軸のいずれかに平行な線を一本引く。この線からの距離がD以下となるような点の個数を最大化したい。最適な線を引いた場合にえられるこの様な点の個数の最大値を求めよ。
　0 ≦ N ≦ 50
　-1000000 ≦ x_i ≦ 1000000
　-1000000 ≦ y_i ≦ 1000000

方針

　まず、x-軸に平行な線を引く場合を考える。点を2つ選んでその2つの距離が2*D以下なら、その2つの点および2つの点の間にある点すべてを距離D以内に含むような線が引ける。これは幾何の問題で良く使う”ギリギリ”を考える手法。
　y-軸も同様に可能。対角軸に平行な線を考えるのはちょっと面倒。点(x1, y1)を通りx+y=0に平行な直線（L1と定義する）と点(x2, y2)を通りx+y=0に平行な直線（L2の定義する）の距離を計算する場合は、L1とx+y=0の距離をまず測り、次にL2とx+y=0の距離を測って、その差分を取ればよいことが分かる。

　あとは、平方根が出てくるので丸め誤差に注意。double使いたくない場合は、長さの2乗を考えればOK。

ソースコード

　md[]とかcd[]とか使ってるけど、これはmain diagonalとcollateral diagonalの略。それぞれ行列の対角成分を呼ぶ呼び方らしい（今日別の問題解いてるときに知った）。main diagonalは右下がりのいわゆる主対角成分。collateral diagonalは右上がりの対角成分。

view plainprint?

1. #define REP(i, n) for(int i=0; i<(int)(n); i++)  
2. #define FOR(i, s, e) for (int i = (int)(s); i < (int)(e); i++)  
3.   
4. class PowerSupply {  
5.     int maxLen(vector<int> x, long long len) {  
6.         int ret = 0;  
7.   
8.         sort(x.begin(), x.end());  
9.         REP(i, x.size()) {  
10.             FOR (j, i, x.size()) {  
11.                 long long diff = x[j] - x[i];  
12.                 diff *= diff;  
13.                 if (diff <= len)  
14.                     ret = max(ret, j-i+1);  
15.             }  
16.         }  
17.   
18.         return ret;  
19.     }  
20.   
21. public:  
22.     int maxProfit(vector<int> x, vector<int> y, int D) {  
23.         int ret = 0;  
24.         int N = x.size();  
25.   
26.         vector<int> md = x;  
27.         vector<int> cd = x;  
28.   
29.         REP(i, N) md[i] = x[i] + y[i];  
30.         REP(i, N) cd[i] = x[i] - y[i];  
31.   
32.         long long D1 = 4LL * D * D;           // orthogonal  
33.         long long D2 = 8LL * D * D;              // diagonal  
34.   
35.         ret = max(ret, maxLen(x, D1));  
36.         ret = max(ret, maxLen(y, D1));  
37.         ret = max(ret, maxLen(md, D2));  
38.         ret = max(ret, maxLen(cd, D2));  
39.   
40.         return ret;  
41.     }  
42. };  

#define REP(i, n) for(int i=0; i<(int)(n); i++)
#define FOR(i, s, e) for (int i = (int)(s); i < (int)(e); i++)

class PowerSupply {
    int maxLen(vector<int> x, long long len) {
        int ret = 0;

        sort(x.begin(), x.end());
        REP(i, x.size()) {
            FOR (j, i, x.size()) {
                long long diff = x[j] - x[i];
                diff *= diff;
                if (diff <= len)
                    ret = max(ret, j-i+1);
            }
        }

        return ret;
    }

public:
    int maxProfit(vector<int> x, vector<int> y, int D) {
        int ret = 0;
        int N = x.size();

        vector<int> md = x;
        vector<int> cd = x;

        REP(i, N) md[i] = x[i] + y[i];
        REP(i, N) cd[i] = x[i] - y[i];

        long long D1 = 4LL * D * D;           // orthogonal
        long long D2 = 8LL * D * D;              // diagonal

        ret = max(ret, maxLen(x, D1));
        ret = max(ret, maxLen(y, D1));
        ret = max(ret, maxLen(md, D2));
        ret = max(ret, maxLen(cd, D2));

        return ret;
    }
};

Google Code Jam Japan 2011 決勝 バクテリアの増殖

解説読んでも分からなかった問題でしたが、再チャレンジしてみたら解けました。この問題は去年一番印象に残った問題の一つです。

問題概要

　あるバクテリアは、現在x個あったとすると１時間後にはx^x個に増えるという性質をもつ。もともとB個いたバクテリアがA時間後には何個になっているかを求めよ。ただし、とても大きな数なのでCで割った余りを求めよ。

1 ≦ A, B, C ≦ 1000

方針

　べき乗の部分は、先に余りを求めて計算してはいけないことに注意。Cが高々1000なので、鳩ノ巣原理によって任意の正の整数pのべき乗を並べたものは最高でも1000回以内にループに入る。ループに入るまでの長さと、ループの１サイクルの周期を予め計算しておく。あとはこの情報を利用して再帰的に解く。実装は結構ややこしかった。printfデバッグがうまくなった気がする。

ソースコード

view plainprint?

1. const int MAX_N = 1000;  
2. const long long INF = 10000;  
3.   
4. int head[MAX_N+1][MAX_N+1];  
5. int cycle[MAX_N+1][MAX_N+1];  
6. int dp[MAX_N+1][MAX_N+1];  
7. int val[MAX_N+1];  
8.   
9. void init() {  
10.     int used[MAX_N+1];  
11.     for (int i = 1; i <= MAX_N; i++) {  
12.         for (int j = 1; j <= MAX_N; j++) {  
13.             memset(used, -1, sizeof(used));  
14.             int p = 0;  
15.             int x = 1;  
16.             while (1) {  
17.                 x = x * i % j;  
18.                 if (used[x] >= 0) {  
19.                     head[i][j] = used[x];  
20.                     cycle[i][j] = p - used[x];  
21.                     break;  
22.                 }  
23.                 used[x] = p++;  
24.             }  
25.         }  
26.     }  
27.     cerr << "init() done!" << endl;  
28. }  
29.   
30. long long power(long long x, long long b, long long mod) {  
31.     long long ret = 1LL;  
32.   
33.     while (b > 0) {  
34.         if (b & 1)  
35.             ret = ret * x % mod;  
36.         x = x * x % mod;  
37.         b >>= 1;  
38.     }  
39.   
40.     return ret;  
41. }  
42.   
43. long long power(long long x, long long b) {  
44.     long long ret = 1LL;  
45.   
46.     while (b > 0) {  
47.         if (b & 1)  
48.             ret = min(INF, ret*x);  
49.         x = min(INF, x*x);  
50.         b >>= 1;  
51.     }  
52.   
53.     return ret;  
54. }  
55.   
56. int rec(int x, int b, int mod) {  
57.     if (!b) return x % mod;  
58.     if (mod == 1) return 0;  
59.     if (dp[b][mod] >= 0) return dp[b][mod];  
60.   
61.     int y = rec(x, b-1, mod);  
62.     if (!y)  
63.         return dp[b][mod] = 0;  
64.   
65.     int p = head[y][mod];  
66.     int q = cycle[y][mod];  
67.   
68.     int ret;  
69.     if (val[b-1] <= p)  
70.         ret = power(y, val[b-1], mod);  
71.     else {  
72.         int z = rec(x, b-1, q);  
73.         z = ((z - (p+1)) % q + q) % q + 1;  
74.         ret = power(y, p+z, mod);  
75.     }  
76.   
77.     return dp[b][mod] = ret;  
78. }  
79.   
80. void solve() {  
81.     int T;  
82.   
83.     cin >> T;  
84.     for (int i = 0; i < T; i++) {  
85.         cerr << "solving case " << i+1 << endl;  
86.         int A, B, C;  
87.         cin >> A >> B >> C;  
88.   
89.         val[0] = A;  
90.         for (int j = 1; j <= MAX_N; j++)  
91.             val[j] = power(val[j-1], val[j-1]);  
92.   
93.         memset(dp, -1, sizeof(dp));  
94.         int ret = rec(A, B, C);  
95.         cout << "Case #" << i+1 << ": " << ret << endl;  
96.     }  
97. }  
98.   
99. #define SUBMIT  
100. int main() {  
101. #ifdef SUBMIT  
102.     freopen("./test.in", "r", stdin);  
103.     freopen("./test.out", "w", stdout);  
104. #endif  
105.   
106.     init();  
107.     solve();  
108.   
109. #ifdef SUBMIT  
110.     fflush(stdout);  
111.     cerr << "all done!" << endl;  
112. #endif  
113.     return 0;  
114. }  

const int MAX_N = 1000;
const long long INF = 10000;

int head[MAX_N+1][MAX_N+1];
int cycle[MAX_N+1][MAX_N+1];
int dp[MAX_N+1][MAX_N+1];
int val[MAX_N+1];

void init() {
    int used[MAX_N+1];
    for (int i = 1; i <= MAX_N; i++) {
        for (int j = 1; j <= MAX_N; j++) {
            memset(used, -1, sizeof(used));
            int p = 0;
            int x = 1;
            while (1) {
                x = x * i % j;
                if (used[x] >= 0) {
                    head[i][j] = used[x];
                    cycle[i][j] = p - used[x];
                    break;
                }
                used[x] = p++;
            }
        }
    }
    cerr << "init() done!" << endl;
}

long long power(long long x, long long b, long long mod) {
    long long ret = 1LL;

    while (b > 0) {
        if (b & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        b >>= 1;
    }

    return ret;
}

long long power(long long x, long long b) {
    long long ret = 1LL;

    while (b > 0) {
        if (b & 1)
            ret = min(INF, ret*x);
        x = min(INF, x*x);
        b >>= 1;
    }

    return ret;
}

int rec(int x, int b, int mod) {
    if (!b) return x % mod;
    if (mod == 1) return 0;
    if (dp[b][mod] >= 0) return dp[b][mod];

    int y = rec(x, b-1, mod);
    if (!y)
        return dp[b][mod] = 0;

    int p = head[y][mod];
    int q = cycle[y][mod];

    int ret;
    if (val[b-1] <= p)
        ret = power(y, val[b-1], mod);
    else {
        int z = rec(x, b-1, q);
        z = ((z - (p+1)) % q + q) % q + 1;
        ret = power(y, p+z, mod);
    }

    return dp[b][mod] = ret;
}

void solve() {
    int T;

    cin >> T;
    for (int i = 0; i < T; i++) {
        cerr << "solving case " << i+1 << endl;
        int A, B, C;
        cin >> A >> B >> C;

        val[0] = A;
        for (int j = 1; j <= MAX_N; j++)
            val[j] = power(val[j-1], val[j-1]);

        memset(dp, -1, sizeof(dp));
        int ret = rec(A, B, C);
        cout << "Case #" << i+1 << ": " << ret << endl;
    }
}

#define SUBMIT
int main() {
#ifdef SUBMIT
    freopen("./test.in", "r", stdin);
    freopen("./test.out", "w", stdout);
#endif

    init();
    solve();

#ifdef SUBMIT
    fflush(stdout);
    cerr << "all done!" << endl;
#endif
    return 0;
}