TopCoder SRM 401 DIV1 Hard - NCool (本番 5 人)

最近、ABC 201〜300 の D 問題埋めを推奨している身としては、僕も同様のトレーニングとして SRM 401〜650 辺りの DIV1 Hard 埋めを始めてみようと思い立った。

• 問題へのリンク
• editorial へのリンク

問題概要

二次元平面上において、以下の条件を満たす線分の cool 度は であるという。

• 両端点が格子点である
• 線分上の両端および内部にちょうど 個の格子点が含まれる

ここで、凸多角形が与えられる。凸多角形に内包される、cool 度が 以上である線分をすべて考えて、その線分の両端をなす格子点に色を塗っていく。色を塗られる格子点の個数を求めよ。

(この図は の場合)

制約

• 多角形の頂点数は 以上 以下
• 多角形の各頂点の座標は 以上 以下
• 

考えたこと

この問題では cool 度が 以上である線分をすべて考えてよいとあるが、実際には cool 度がちょうど である線分のみを考えれば十分。その方が数えやすそうだ。

ここで、cool 度が であるという条件を言い換えていくことを考える。これはもう有名問題で結論が出ていて、線分の cool 度が であるとは、

「線分の両端点の x 座標の差と y 座標の差の最大公約数が 」

となる。ということは、cool 度が である線分の両端点を とすると

• 
• 

が成り立つ。これは必要条件であって十分条件でないが、とても扱いやすそうな条件なので、この条件が成り立つような線分の両端点に色を塗っていけば十分だと言えたら嬉しい！！！

, が成り立つとき、線分の cool 度は ( は正の整数) と表せることは言える。もともと問題は「cool 度が 以上である線分」を考えていたくらいなので、「cool 度が と表せる線分」をすべて考えることにしても全然 OK！！！

帰着した問題

以上の考察から、次の問題へと帰着された。

---

各格子点 を、 に基づいて、 通りに分類する。

与えられた凸多角形内部に、同じグループに分類される格子点が複数あるならば、それらの格子点に色を塗ることにする。

色が塗られた格子点の個数を求めよ。

---

ここまで来れば、あとはなんとでもなると思われた。まず、x 座標や y 座標の値としてありうる値の最大値を (= 10000) とする。 の場合は 0 と返してよい。

よって、これ以降は と仮定できる。

座標を で割った余り を固定して考える。 に対して、その 座標をもち、凸多角形の内部にあるような格子点の 座標の最小値と最大値を求める方針をとった。

適宜、いもす法なども活用することで、次の配列が の計算量で求められることがわかった。

---

num[i][j] ← 凸多角形に内包される、 座標と 座標を で割った余りがそれぞれ であるような格子点の個数

---

全体の計算量は と評価できる。

コード

moj プラグインに基づくテストも含む。

#include <bits/stdc++.h>
using namespace std;
using pint = pair<int, int>;
using pll = pair<long long, long long>;
template<class T> inline bool chmax(T& a, T b) { if (a < b) { a = b; return 1; } return 0; }
template<class T> inline bool chmin(T& a, T b) { if (a > b) { a = b; return 1; } return 0; }

#define REP(i, n) for (long long i = 0; i < (long long)(n); ++i)
#define REP2(i, a, b) for (long long i = a; i < (long long)(b); ++i)
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl
template<class T1, class T2> ostream& operator << (ostream &s, pair<T1,T2> P)
{ return s << '<' << P.first << ", " << P.second << '>'; }
template<class T> ostream& operator << (ostream &s, vector<T> P)
{ for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }
template<class T> ostream& operator << (ostream &s, deque<T> P)
{ for (int i = 0; i < P.size(); ++i) { if (i > 0) { s << " "; } s << P[i]; } return s; }
template<class T> ostream& operator << (ostream &s, vector<vector<T> > P)
{ for (int i = 0; i < P.size(); ++i) { s << endl << P[i]; } return s << endl; }
template<class T> ostream& operator << (ostream &s, set<T> P)
{ for(auto it : P) { s << "<" << it << "> "; } return s; }
template<class T> ostream& operator << (ostream &s, multiset<T> P)
{ for(auto it : P) { s << "<" << it << "> "; } return s; }
template<class T1, class T2> ostream& operator << (ostream &s, map<T1,T2> P)
{ for(auto it : P) { s << "<" << it.first << "->" << it.second << "> "; } return s; }


/*/////////////////////////////*/
// 幾何ライブラリ
/*/////////////////////////////*/

// basic settings
using DD = long double;
constexpr long double PI = 3.141592653589793238462643383279502884L;
constexpr long double INF = 1LL<<60;  // to be set appropriately
constexpr long double EPS = 1e-10;    // to be set appropriately
long double torad(int deg) {return (long double)(deg) * PI / 180;}
long double todeg(long double ang) {return ang * 180 / PI;}

// Point or Vector
struct Point {
    DD x, y;
    
    // constructor
    constexpr Point() : x(0), y(0) {}
    constexpr Point(DD x, DD y) : x(x), y(y) {}
    
    // various functions
    constexpr Point conj() const {return Point(x, -y);}
    constexpr DD dot(const Point &r) const {return x * r.x + y * r.y;}
    constexpr DD cross(const Point &r) const {return x * r.y - y * r.x;}
    constexpr DD norm() const {return dot(*this);}
    constexpr long double abs() const {return sqrt(norm());}
    constexpr long double amp() const {
        long double res = atan2(y, x);
        if (res < 0) res += PI*2;
        return res;
    }
    constexpr bool eq(const Point &r) const {return (*this - r).abs() <= EPS;}
    constexpr Point rot90() const {return Point(-y, x);}
    constexpr Point rot(long double ang) const {
        return Point(cos(ang) * x - sin(ang) * y, sin(ang) * x + cos(ang) * y);
    }
    
    // arithmetic operators
    constexpr Point operator - () const {return Point(-x, -y);}
    constexpr Point operator + (const Point &r) const {return Point(*this) += r;}
    constexpr Point operator - (const Point &r) const {return Point(*this) -= r;}
    constexpr Point operator * (const Point &r) const {return Point(*this) *= r;}
    constexpr Point operator / (const Point &r) const {return Point(*this) /= r;}
    constexpr Point operator * (DD r) const {return Point(*this) *= r;}
    constexpr Point operator / (DD r) const {return Point(*this) /= r;}
    constexpr Point& operator += (const Point &r) {
        x += r.x, y += r.y;
        return *this;
    }
    constexpr Point& operator -= (const Point &r) {
        x -= r.x, y -= r.y;
        return *this;
    }
    constexpr Point& operator *= (const Point &r) {
        DD tx = x, ty = y;
        x = tx * r.x - ty * r.y;
        y = tx * r.y + ty * r.x;
        return *this;
    }
    constexpr Point& operator /= (const Point &r) {
        return *this *= r.conj() / r.norm();
    }
    constexpr Point& operator *= (DD r) {
        x *= r, y *= r;
        return *this;
    }
    constexpr Point& operator /= (DD r) {
        x /= r, y /= r;
        return *this;
    }

    // friend functions
    friend ostream& operator << (ostream &s, const Point &p) {
        return s << '(' << p.x << ", " << p.y << ')';
    }
    friend constexpr Point conj(const Point &p) {return p.conj();}
    friend constexpr DD dot(const Point &p, const Point &q) {return p.dot(q);}
    friend constexpr DD cross(const Point &p, const Point &q) {return p.cross(q);}
    friend constexpr DD norm(const Point &p) {return p.norm();}
    friend constexpr long double abs(const Point &p) {return p.abs();}
    friend constexpr long double amp(const Point &p) {return p.amp();}
    friend constexpr bool eq(const Point &p, const Point &q) {return p.eq(q);}
    friend constexpr Point rot90(const Point &p) {return p.rot90();}
    friend constexpr Point rot(const Point &p, long long ang) {return p.rot(ang);}
};

// necessary for some functions
constexpr bool operator < (const Point &p, const Point &q) {
    return (abs(p.x - q.x) > EPS ? p.x < q.x : p.y < q.y);
}

// Line
struct Line : vector<Point> {
    Line(Point a = Point(0.0, 0.0), Point b = Point(0.0, 0.0)) {
        this->push_back(a);
        this->push_back(b);
    }
    friend ostream& operator << (ostream &s, const Line &l) {
        return s << '{' << l[0] << ", " << l[1] << '}';
    }
};

int ccw_for_dis(const Point &a, const Point &b, const Point &c) {
    if (cross(b-a, c-a) > EPS) return 1;
    if (cross(b-a, c-a) < -EPS) return -1;
    if (dot(b-a, c-a) < -EPS) return 2;
    if (norm(b-a) < norm(c-a) - EPS) return -2;
    return 0;
}
Point proj(const Point &p, const Line &l) {
    DD t = dot(p - l[0], l[1] - l[0]) / norm(l[1] - l[0]);
    return l[0] + (l[1] - l[0]) * t;
}
Point refl(const Point &p, const Line &l) {
    return p + (proj(p, l) - p) * 2;
}
bool is_inter_PL(const Point &p, const Line &l) {
    return (abs(p - proj(p, l)) < EPS);
}
bool is_inter_PS(const Point &p, const Line &s) {
    return (ccw_for_dis(s[0], s[1], p) == 0);
}
bool is_inter_LL(const Line &l, const Line &m) {
    return (abs(cross(l[1] - l[0], m[1] - m[0])) > EPS ||
            abs(cross(l[1] - l[0], m[0] - l[0])) < EPS);
}
bool is_inter_SS(const Line &s, const Line &t) {
    if (eq(s[0], s[1])) return is_inter_PS(s[0], t);
    if (eq(t[0], t[1])) return is_inter_PS(t[0], s);
    return (ccw_for_dis(s[0], s[1], t[0]) * ccw_for_dis(s[0], s[1], t[1]) <= 0 &&
            ccw_for_dis(t[0], t[1], s[0]) * ccw_for_dis(t[0], t[1], s[1]) <= 0);
}
DD distance_PL(const Point &p, const Line &l) {
    return abs(p - proj(p, l));
}
DD distance_PS(const Point &p, const Line &s) {
    Point h = proj(p, s);
    if (is_inter_PS(h, s)) return abs(p - h);
    return min(abs(p - s[0]), abs(p - s[1]));
}
DD distance_LL(const Line &l, const Line &m) {
    if (is_inter_LL(l, m)) return 0;
    else return distance_PL(m[0], l);
}
DD distance_SS(const Line &s, const Line &t) {
    if (is_inter_SS(s, t)) return 0;
    else return min(min(distance_PS(s[0], t), distance_PS(s[1], t)),
                    min(distance_PS(t[0], s), distance_PS(t[1], s)));
}

Point proj_for_crosspoint(const Point &p, const Line &l) {
    DD t = dot(p - l[0], l[1] - l[0]) / norm(l[1] - l[0]);
    return l[0] + (l[1] - l[0]) * t;
}
vector<Point> crosspoint(const Line &l, const Line &m) {
    vector<Point> res;
    DD d = cross(m[1] - m[0], l[1] - l[0]);
    if (abs(d) < EPS) return vector<Point>();
    res.push_back(l[0] + (l[1] - l[0]) * cross(m[1] - m[0], m[1] - l[0]) / d);
    return res;
}
vector<Point> crosspoint_SS(const Line &l, const Line &m) {
    if (is_inter_SS(l, m)) return crosspoint(l, m);
    else return vector<Point>();
}

// 凸包 (一直線上の3点を含めない)
vector<Point> convex_hull(vector<Point> &ps) {
    int n = (int)ps.size();
    vector<Point> res(2*n);
    auto cmp = [&](Point p, Point q) -> bool {
        return (abs(p.x - q.x) > EPS ? p.x < q.x : p.y < q.y);
    };
    sort(ps.begin(), ps.end(), cmp);
    int k = 0;
    for (int i = 0; i < n; ++i) {
        if (k >= 2) {
            while (cross(res[k-1] - res[k-2], ps[i] - res[k-2]) < EPS) {
                --k;
                if (k < 2) break;
            }
        }
        res[k] = ps[i]; ++k;
    }
    int t = k+1;
    for (int i = n-2; i >= 0; --i) {
        if (k >= t) {
            while (cross(res[k-1] - res[k-2], ps[i] - res[k-2]) < EPS) {
                --k;
                if (k < t) break;
            }
        }
        res[k] = ps[i]; ++k;
    }
    res.resize(k-1);
    return res;
}
/* 幾何ライブラリここまで */

const int MAX = 10000;
// 座標 x において、多角形の外周または内部に含まれる点の y 座標の最小値と最大値を求める
pair<long long, long long> calc(const vector<Point> &vp, int x) {
    DD minv = INF, maxv = -INF;
    for (int i = 0; i < vp.size(); ++i) {
        Line seg(vp[i], vp[(i+1)%vp.size()]);
        Line tate(Point(x, -1), Point(x, MAX + 1));
        const auto &inter = crosspoint_SS(seg, tate);
        if (inter.empty()) continue;
        chmin(minv, inter[0].y);
        chmax(maxv, inter[0].y);
    }
    long long mi = (long long)(minv - EPS + 1);
    long long ma = (long long)(maxv + EPS);
    return {mi, ma};
}


class NCool {
public:
    int nCoolPoints(vector <int> x, vector <int> y, int N) {
        --N;
        if (N > MAX) return 0;
        
        vector<Point> vp(x.size());
        for (int i = 0; i < x.size(); ++i) vp[i] = Point(x[i], y[i]);
        vp = convex_hull(vp);
        
        long long res = 0;
        for (int i = 0; i < N; ++i) {
            vector<long long> num(N+1, 0);
            
            //cout << "-----------------=" << endl; COUT(i);
            
            // 0 以上 r 未満の値について、num に d 足す
            auto add = [&](int r, long long d) -> void {
                long long q = r / N;
                r %= N;
                num[0] += q * d, num[N] -= q * d;
                num[0] += d, num[r] -= d;
            };
            
            for (int x = i; x <= MAX; x += N) {
                auto [mi, ma] = calc(vp, x);
                if (ma < 0) continue;
                
                //cout << x << ": " << mi << ", " << ma << endl;
                
                add(ma+1, 1), add(mi, -1);
            }
            for (int i = 0; i < N; ++i) num[i+1] += num[i];
            
            //COUT(num);
            for (int i = 0; i < N; ++i) if (num[i] > 1) res += num[i];
        }
        return res;
    }
};



// BEGIN CUT HERE
namespace moj_harness {
    int run_test_case(int);
    void run_test(int casenum = -1, bool quiet = false) {
        if (casenum != -1) {
            if (run_test_case(casenum) == -1 && !quiet) {
                cerr << "Illegal input! Test case " << casenum << " does not exist." << endl;
            }
            return;
        }
        
        int correct = 0, total = 0;
        for (int i=0;i <= 10; ++i) {
            int x = run_test_case(i);
            if (x == -1) {
                if (i >= 100) break;
                continue;
            }
            correct += x;
            ++total;
        }
        
        if (total == 0) {
            cerr << "No test cases run." << endl;
        } else if (correct < total) {
            cerr << "Some cases FAILED (passed " << correct << " of " << total << ")." << endl;
        } else {
            cerr << "All " << total << " tests passed!" << endl;
        }
    }
    
    int verify_case(int casenum, const int &expected, const int &received, clock_t elapsed) { 
        cerr << "Example " << casenum << "... "; 
        
        string verdict;
        vector<string> info;
        char buf[100];
        
        if (elapsed > CLOCKS_PER_SEC / 200) {
            sprintf(buf, "time %.2fs", elapsed * (1.0/CLOCKS_PER_SEC));
            info.push_back(buf);
        }
        
        if (expected == received) {
            verdict = "PASSED";
        } else {
            verdict = "FAILED";
        }
        
        cerr << verdict;
        if (!info.empty()) {
            cerr << " (";
            for (int i=0; i<(int)info.size(); ++i) {
                if (i > 0) cerr << ", ";
                cerr << info[i];
            }
            cerr << ")";
        }
        cerr << endl;
        
        if (verdict == "FAILED") {
            cerr << "    Expected: " << expected << endl; 
            cerr << "    Received: " << received << endl; 
        }
        
        return verdict == "PASSED";
    }

    int run_test_case(int casenum__) {
        switch (casenum__) {
            case 0: {
                int x[]                   = {0, 1, 2, 7, 7};
                int y[]                   = {3, 1, 6, 1, 5};
                int n                     = 6;
                int expected__            = 21;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 1: {
                int x[]                   = {0, 1, 0};
                int y[]                   = {0, 0, 1};
                int n                     = 2;
                int expected__            = 3;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 2: {
                int x[]                   = {0, 0, 1, 2, 2, 1, 0, 0, 2};
                int y[]                   = {0, 1, 2, 2, 1, 0, 0, 0, 2};
                int n                     = 3;
                int expected__            = 6;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 3: {
                int x[]                   = {0, 1, 1, 2, 2, 3, 3, 4, 4, 5};
                int y[]                   = {1, 0, 2, 0, 2, 0, 2, 0, 2, 1};
                int n                     = 5;
                int expected__            = 4;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 4: {
                int x[]                   = {0, 1, 1, 2, 2, 3, 3, 4, 4, 5};
                int y[]                   = {1, 0, 2, 0, 2, 0, 2, 0, 2, 1};
                int n                     = 4;
                int expected__            = 10;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
                
                // custom cases
                
            case 5: {
                int x[]                   = {0, 10000, 10000, 0};
                int y[]                   = {0, 0, 10000, 10000};
                int n                     = 10001;
                int expected__            = 40000;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 6: {
                int x[]                   = {0, 10000, 0};
                int y[]                   = {0, 0, 10000};
                int n                     = 10001;
                int expected__            = 3;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            case 7: {
                int x[]                   = {0, 10000, 10000, 0};
                int y[]                   = {0, 0, 10000, 10000};
                int n                     = 10002;
                int expected__            = 0;
                
                clock_t start__           = clock();
                int received__            = NCool().nCoolPoints(vector <int>(x, x + (sizeof x / sizeof x[0])), vector <int>(y, y + (sizeof y / sizeof y[0])), n);
                return verify_case(casenum__, expected__, received__, clock()-start__);
            }
            default:
                return -1;
        }
    }
}
 

int main(int argc, char *argv[]) {
    if (argc == 1) {
        moj_harness::run_test();
    } else {
        for (int i=1; i<argc; ++i)
            moj_harness::run_test(atoi(argv[i]));
    }
}
// END CUT HERE

TopCoder SRM 400 DIV1 Hard - CollectingBonuses (本番 2 人)

ついこの間の AGC でもあった「〜が終了するまでの回数の期待値」に関する問題！！！

それにしても数値誤差しゃん...
150 人が提出してAC 2 人という恐ろしい問題。
　
実はこの回の writer であった ymatux さんから、直接話を聞いたのだった。Petr はこの回 resubmit をしていて、148 人が WA った中 AC を勝ち取ったらしい。「やっぱり Petr は違うなーと思った」と言っていた。

問題概要

種類の品目が当確率で排出されるガチャがある。このガチャで 種類の異なる品目を揃えるまでに引く回数の期待値を、相対誤差または絶対誤差が 以内の範囲で求めよ。

制約

• 

考えたこと

答えの数式を得ることはそれほど難しくない

• 1 個目を得るまでの期待値: 1 回
• 2 個目を得るまでの期待値: 回
• 3 個目を得るまでの期待値: 回
• ...
• 個目を得るまでの期待値: 回

というわけなのでこれを合計すればよい。もっというと

• 

として

• 

が答えである。難しいのは「」という制約と「数値誤差」。数値誤差の怖さに気づかずに大量の Failed System Error が出たようだ。

H(n) の計算: n <= 10¹⁸ なので

がでかいので愚直に全部足すことはできない。なので解析的になんとかしよう。

はいわゆる調和級数とよばれるものなので、解析的な取り扱いができそうである。数IIIで習った区分求積法によって、以下のことがいえる。

• 

ここで はオイラー定数とよばれる値で、0.57... とかなんとか。これの収束性が気になるが、以下の記述がある！

これによるとなんと、

• 

とした方が近似精度がよい！！！そして二乗の精度だとわかる。なので の精度で求めたかったら、

• くらいまでは自力で計算
• では を計算

という感じで を求めることができる

まだ罠があるっぽい！！！

ところがこれでもまだ罠があって「そんなん気づくかーーー！！！」という感じなのだけれども、 　

　

が答えだけど、 を安易に計算してはいけない！！！一般に数値計算において値が近い浮動小数同士の引き算は超危険で、たとえば 　

1.3244324434... - 1.3244324431... = 0.00000000003... 　

となって、相対精度がどっかーんと落ちてしまう。これを回避するために、 が十分小さいときは 　

　

とすることが考えられる。 がある程度大きければこう考えて良いと思われるので作戦としては

• の計算の先頭部分について、分母が十分大きくなるまで (安全のため くらい) やる
• 残りは、 を活用

という感じでやろう。

まだまだ罠があるっぽい！！！

ここまでケアしてもまだ罠があって、 の の中身は に近いので、このままだと精度が低い。そこで、log1p を用いるとよい。これは を計算してくれる関数。 というわけで 　 　

ということになる。

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <cmath>
using namespace std;
#define COUT(x) cout << #x << " = " << (x) << " (L" << __LINE__ << ")" << endl

// n/n + n/(n-1) + n/(n-2) + ... + n/(n-k+1)
// f(n) = 1 + 1/2 + ... + 1/n として
// n(f(n) - f(n-k)) が答え

class CollectingBonuses {
public:
    double expectedBuy(string sn, string sk) {
        long double res = 0.0;
        const long long M = 100000000;

        istringstream isn(sn); long long n; isn >> n;
        istringstream isk(sk); long long k; isk >> k;

        long long m = n - k + 1;
        while (m < M) {
            res += (long double)(1.0) / m;
            if (m == n) return res * n;
            ++m;
        }
        --m;

        // f(n) - f(m) 〜 log((2n+1)/(2m+1)) = log(1 + 2(n-m)/(2m+1))
        res += log1pl((n - m) * 2.0 / (m * 2.0 + 1));
        return res * n;
    }
};



TopCoder SRM 309 DIV1 Hard - StoneGameStrategist (本番 4 人)

staircase nim の流れで

問題へのリンク

問題概要

個の石の山が左から順に一列に並んでいて、各山には 個の石が積まれている。初期状態では

• 

を満たしている。今先手と後手が交互に

• 石が 1 個以上ある好きな山を 1 つ選んで
• 何個かの石を取り去る
• ただし という状態を維持しなければならない、この条件を満たす限りは任意個数の石を取り去ることができる

という条件でプレイする。先に操作を行えなくなった方が負けである。双方最善を尽くしたとき、どちらが勝つか？

また先手必勝の場合は、最初にとるべき石の個数として考えられる最小値を求めよ。

制約

• 

考えたこと

実は

• 隣の山に石を渡していく Nim

と同じ問題になっている。具体的には、例えば

{3, 9, 10} -> {3, 7, 10}

のように石を減らす操作は、隣り合う石の個数の差でみると

{3, 6, 1} -> {3, 4, 3}

というように、「左隣の山から石を何個か選んで右隣の山に石を移す操作」と思うことができる。そして「単調非減少」という条件は、石の個数の差の世界では「どの山の石の個数もマイナスにはならないように」というのにピッタリ対応する。

よって隣り合う石の個数で見た世界において、必勝法は

• ^ ^ = 0

となることであり、実際に取るべき石の個数も比較的容易に復元できる。

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
using namespace std;

class StoneGameStrategist {
public:
    string play(vector<int> a) {
        vector<int> b;
        b.push_back(a[0]);
        for (int i = 0; i + 1 < a.size(); ++i) {
            b.push_back(a[i+1] - a[i]);
        }
        int nim = 0;
        for (int i = b.size() - 1; i >= 0; i -= 2) nim ^= b[i];
        if (nim == 0) return "YOU LOSE";

        int mi = 1<<29, mii = -1;
        // i を試す場合を考える
        for (int i = b.size() - 1; i >= 0; i -= 2) {
            int need = nim ^ b[i];

            // i から i+1 へ移す、これは a では i から取ることに相当
            if (need < b[i]) {
                if (mi >= b[i] - need) mi = b[i] - need, mii = i;
            }
            // i-1 から i へ移す、これは a では i-1 から取ることに相当
            else {
                // 左がマイナスになってはダメ
                if (i == 0 || need - b[i] > b[i-1]) continue;
                if (mi >= need - b[i]) mi = need - b[i], mii = i-1;
            }
        }
        stringstream ss;
        ss << "TAKE " << mi << " STONES FROM PILE " << mii;
        return ss.str();
    }
};

TopCoder SRM 735 DIV1 Medium - QuadraticIdentity

ちょうど中国剰余定理シリーズやっていたので、ピンポイントだった。

問題概要

整数 が与えられる。

(mod. )

を満たすような < ) をすべて求め、そのサイズが 500 より大きい場合には 500 以下になるまで以下の操作をした上で出力せよ:

• 数列を小さい順にソートして奇数番目のみを残す

制約

• 

解法

サイズが 500 より大きかったら...という話が怪しい雰囲気だけど、単に topcoder の出力受け取りキャパシティの問題っぽい？？？全然問題とは関係ない感じっぽい。

さて、

(mod. )

になる。 は の倍数というわけだが、x と x-1 とは互いに素なので、

と素因数分解したときに、 と とは共通の素因子をもたない。したがって、 個の素数べき を と とに振り分ける感じになる。その方法は 通りある。振り分けて

• (mod. )
• (mod. )

になったとする。これを満たす (mod. ) は、中国剰余定理を用いて求めることができる。ライブラリで殴ってもいいのだが、もっと簡単に求めることもできる。

mod. での の逆元を として、

• 

としてあげると、 (mod. ) と、 (mod. ) をともに満たす。

最後に素因数分解したときの素因子の種類数 がどの程度の大きさになり得るかを見積もっておく。Python でサッと計算すると

>>> 2*3*5*7*11*13*17*19*23*29*31*37*41
304250263527210
>>> 2*3*5*7*11*13*17*19*23*29*31*37*41*43
13082761331670030

となったので、最悪でも であることがわかる。したがって、求める の個数は最悪でも 個であり、大分余裕がある。

#include <iostream>
#include <vector>
#include <algorithm>
using namespace std;

vector<long long> prime_factor(long long n) {
    vector<long long> res;
    for (long long i = 2; i*i <= n; ++i) {
        long long num = 1;
        if (n % i == 0) {
            while (n%i == 0) {
                num *= i;
                n /= i;
            }
            res.push_back(num);
        }
    }
    if (n != 1) res.push_back(n);
    return res;
}

inline long long mod(long long a, long long m) {
    return (a % m + m) % m;
}

inline long long inv(long long a, long long m) {
    long long b = m, u = 1, v = 0;
    while (b) {
        long long t = a/b;
        a -= t*b; swap(a, b);
        u -= t*v; swap(u, v);
    }
    return mod(u, m);
}

class QuadraticIdentity {
public:
    vector<long long> getFixedPoints(long long m) {
        vector<long long> pf = prime_factor(m);
        vector<long long> list;
        int n = (int)pf.size();
        for (int bit = 0; bit < (1<<n); ++bit) {
            long long A = 1, B = 1;
            for (int i = 0; i < pf.size(); ++i) {
                if (bit & (1<<i)) A *= pf[i];
                else B *= pf[i];
            }
            long long tmp = inv(A, B) * A;
            list.push_back(tmp);
        }
        sort(list.begin(), list.end());
        
        while (list.size() > 500) {
            vector<long long> nlist;
            for (int i = 0; i < (int)list.size(); i += 2) nlist.push_back(list[i]);
            list = nlist;
        }
        
        return list;
    }
};

TopCoder SRM 694 DIV1 Medium - DistinguishableSetDiv1

難しく考えすぎてしまった。
文字列問題をみてもひるまずに解けるようになりたい。

問題概要

m 文字からなる文字列が n 個与えられる。
{0, 1, 2, …, m-1} の部分集合 bit のうち、
どの 2 つの文字列を選んでも、それらの文字列の bit に対応する部分文字列
(文字列 “string” に対して、bit = {0, 2, 3} に対応する文字列は “sri”)
が互いに等しくならないような部分集合 bit の個数を求めよ。

(制約)
・0 <= m <= 20
・0 <= n <= 1000

解法

部分集合 bit が条件を満たさない
⇔ ある 2 つの文字列 s, t が存在して、s[bit] = t[bit] である。
⇔ ある 2 つの文字列 s, t が存在して、bit は、s と t の極大共通部分文字列の部分文字列である。

そこで、nC2 通りの文字列の組み合わせそれぞれについて、
極大共通部分文字列を求めておいて(同じ文字の index を列挙するのみ)、
それらの部分集合がすべて列挙できればよい。

これは、以下のような dp で実現できる。計算量は、O(mn^2 + m2^m)。

dp[ S ] := 部分集合 S が条件を満たさないとき 1、満たすとき 0
として、

dp[S - {i}] |= dp[S]

コード

#include <iostream>
#include <sstream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstring>
#include <string>
#include <vector>
#include <stack>
#include <queue>
#include <deque>
#include <map>
#include <set>
#include <bitset>
#include <numeric>
#include <utility>
#include <iomanip>
#include <algorithm>
#include <functional>
using namespace std;

bool dp[(1<<21)];

class DistinguishableSetDiv1 {
public:
    int count(vector <string> answer) {
        int n = answer.size();
        int m = answer[0].size();
        int res = 0;
        memset(dp, 0, sizeof(dp));
        for (int i = 0; i < n; ++i) for (int j = i+1; j < n; ++j) {
            int bit = 0;
            for (int k = 0; k < m; ++k) {
                if (answer[i][k] == answer[j][k])
                    bit += (1 << k);
            }
            dp[bit] = true;
        }
        for (int bit = (1 << m) - 1; bit >= 0; --bit) {
            if (dp[bit]) {
                for (int t = bit; t; t &= t-1) {
                    int sub = bit ^ (t & -t);
                    dp[sub] = true;
                }
            }
            else ++res;
        }
        return res;
    }
};