ferinの競プロ帳

競プロについてのメモ

2015 ICL, Finals, Div. 1 J. Ceizenpok’s formula (nCk mod m の求め方)

問題ページ
Problem - J - Codeforces

概要

nCk mod m を求める。 1 <= n <= 10^18, 0 <= k <= n, 2 <= m <= 1000000

考えたこと

JAG夏合宿で見た問題
パット見nCk mod m求めるいつものやつだ!と思って制約を見たら真顔になった
n, kが大きすぎるのでn, kまでの階乗を計算しておくことができないしmが素数とは限らないので逆元が求まらないときもあるしわからなくて終了した
解説でlucasの定理とか中国剰余定理とか聞いたけどそのときはチンプンカンプンだった

{ }

中国剰余定理

2元の連立合同式について考える。
\begin{align} x \equiv a_1\ ({\rm mod}\ m_1) \cdots (1) \\ x \equiv a_2\ ({\rm mod}\ m_2) \cdots (2) \end{align} { m_1, m_2 }が互いに素であるとき { 0 \leq x \leq m_1m_2} の範囲に条件を満たす{ x }がただ一つ存在する。このとき、解xは拡張ユークリッドの互除法を利用して求めることができる。式(1)より
\begin{align} x = a_1 + zm_1 \end{align} となる。
\begin{align} x = a_2 = a_1 + zm_1\ ({\rm mod}\ m_2) \end{align} したがって、上式を満たすような zを見つければよい。 \begin{align} z = (a_2-a_1)k_1 (ただし k_1 \equiv m_1^{-1} {\rm mod}\ m_2) \end{align} と置いたとき上式が成り立つことを示す。 \begin{align} x = a_1 + (a_2-a_1)k_1m_1 \\ k_1m_1 \equiv 1 ({\rm mod}\ m_2) より x \equiv a_2 ({\rm mod}\ m_2) \end{align}

{m_1,m_2}は互いに素なことから逆元は存在する。逆元は拡張ユークリッドの式より求めることができるから{k_1}は求められる。したがって以下の式 \begin{align} x = a_1 + (a_2-a_1)k_1m_1 \end{align} より連立合同式の解{x}を求めることができる。

n個の合同式から成る連立合同式についても同様に解{x}を求めることができる。 \begin{align} x = a_i\ ({\rm mod}\ m_i) (1 \leq i \leq n) \end{align} まず{i=1,2}である2式について条件を満たす解を上の方法を用いて求める。すると、 \begin{align} x = a_{1,2}\ ({\rm mod}\ m_1m_2) \end{align} が求まる。次にこの式と{i=3}の式についての連立合同式の解を求める。これを繰り返していくことで最終的に条件を満たす解{x} \begin{align} x = a'\ ({\rm mod}\ m_1m_2\cdots m_n) \end{align} を求めることができる。

modの素因数分解

{C(n, k)\ ({\rm mod} m)}{m}素因数分解する。 \begin{align} m\ =\ m_1m_2\cdots m_r\ =\ p_1^{q_1}p_2^{q_2}\cdots p_r^{q_r} \end{align} このとき{C(n, k)\ ({\rm mod}\ m_i) (1 \leq i \leq r)}を求めることができたとする。するとr個の連立合同式が導出でき、中国剰余定理を用いることにより \begin{align} nCk\ ({\rm mod}\ m_1m_2\cdots m_r = m) \end{align} を求めることができる。

lucasの定理

では、{C(n, k)\ ({\rm mod}\ m_i) (1 \leq i \leq r)}はどのようにして求めればいいのか。これにはlucasの定理を用いる。素数{p}と非負整数{m,n}について
\begin{align} C(n,k) = \prod_{i=0}^{l} C(n_i, k_i)\ ({\rm mod}\ p) \end{align}
が成り立つ。ここで{n_i,k_i}
\begin{align} n = n_{0} + n_{1}p + \cdots + n_{l-1}p^{l-1} \\ k = k_{0} + k_{1}p + \cdots + k_{l-1}p^{l-1} \end{align}
を表す。したがって、modが素数のときは{O({\rm log}n)}で二項係数{C(n,k)}を求めることができる。

lucasの定理の拡張

modが素数の累乗{p^q}のときに二項係数{C(n,k)}を求める。
{r=n-k}とおく。{n_i,k_i,r_i}を上と同様にp進数に展開したときのi桁目と定義する。また、{N_j,K_j,R_j}{j}桁目から{j+q-1}桁目までの部分列と定義する。 \begin{align} N_j = n_{j} + n_{j+1}p + \cdots + n_{j+q-1}p^{q-1} \\ N_j = [n/p^{j}]\ ({\rm mod}\ p^{q}) \end{align} また、{e_{j}}を \begin{align} e_{j} = ([n/p^{j+1}] + [n/p^{j+2} + \cdots])-([r/p^{j+1}] + [r/p^{j+2} + \cdots]) \\ -([k/p^{j+1}] + [k/p^{j+2} + \cdots]) \end{align} と定義する。さらに{(k!)_{p}}をk以下でpの倍数でないものの積と定義する。これらを用いると二項係数{C(n,k)}は \begin{align} \frac{C(n,k)}{p^{e_0}} = (\pm 1)^{e_{q-1}} (\frac{(N_{0}!)_p}{(M_{0}!)_p(R_{0}!)_p})(\frac{(N_{1}!)_p}{(M_{1}!)_p(R_{1}!)_p}) \cdots (\frac{(N_{l}!)_p}{(M_{l}!)_p(R_{l}!)_p})\ ({\rm mod}\ p^{q}) \end{align} と書ける。{\pm 1}の部分は{p=2, q \geq 3}のときには1、それ以外のときには-1となる。

{e_{j}}{O({\rm log}n + {\rm log}r + {\rm log}k)}で求めることができる。{{k!}_p}はp進数q桁の数字であることを考えると最大でも{p^q}で抑えられる。{O(p^{q})}でとりうる{(k!)_{p}}を列挙する。また、{(k!)_{p}}にはpの倍数が含まれないことからmodの{p^q}とは互いに素である。したがって逆元が存在する。拡張ユークリッドの互除法を用いて逆元を計算することができる。

よってmodが素数の累乗のときの二項係数は計算することができる。

まとめ

まずmodを因数分解し、得られた因数をmodとする二項係数の値をlucasの定理の拡張を用いて求める。その後中国剰余定理で最終的な値を求める。

参考文献

nCr mod mの求め方 - uwicoder - アットウィキ
Coding-Templates/BinCoeff.pdf at master · rishirajsinghjhelumi/Coding-Templates · GitHub
競技プログラミングにおける数学的問題まとめ - はまやんはまやんはまやん
Lucasの定理とその証明 | 高校数学の美しい物語
中国剰余定理の証明と例題(二元の場合) | 高校数学の美しい物語

#define __USE_MINGW_ANSI_STDIO 0
#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
#define int ll
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define PB push_back

const ll LLINF = (1LL<<60);
const int INF = (1LL<<30);
const int MOD = 1000000007;

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }
template<class S,class T>
ostream &operator <<(ostream& out,const pair<S,T>& a){
  out<<'('<<a.first<<','<<a.second<<')';
  return out;
}

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

// ax + by = gcd(a, b) となる {x, y, gcd(a, b)} を返す
// O(log(min(a, b)))
ll extgcd(ll a, ll b, ll &x, ll &y) {
  ll g = a; x = 1, y = 0;
  if(b != 0) g = extgcd(b, a%b, y, x), y -= (a/b) * x;
  return g;
}

// a^-1 mod n を返す 存在しなければ-1
// O(log(n))
ll inv(ll a, ll n) {
  ll s, t;
  extgcd(a, n, s, t);
  return (n+s) % n;
}

// 二分累乗法 x^e % mod O(log(e))
ll binpow(ll x, ll e, ll mod) {
  ll a = 1, p = x;
  while(e > 0) {
    if(e%2 == 0) {p = (p*p) % mod; e /= 2;}
    else {a = (a*p) % mod; e--;}
  }
  return a % mod;
}

// x = a1 mod m1, x2 = a2 mod m2 を解く オーバーフローには注意
// O(log(min(m1, m2)))
pair<ll, ll> crt(ll a1, ll a2, ll m1, ll m2) {
  auto normal = [](ll x, ll m) { return x>=-x ? x%m : m-(-x)%m; };
  auto modmul = [&normal](ll a, ll b, ll m) { return normal(a, m)*normal(b, m)%m; };
  ll k1, k2;
  ll g = extgcd(m1, m2, k1, k2);
  if(normal(a1, g) != normal(a2, g)) return {-1, -1};
  ll l = m1 / g * m2;
  ll x = a1 + modmul(modmul((a2-a1)/g, k1, l), m1, l);
  return {x, l};
}

pair<ll, ll> crt(vector<ll> a, vector<ll> m) {
  ll mod = 1, ans = 0;
  int n = a.size();
  REP(i, n) {
    tie(ans, mod) = crt(ans, a[i], mod, m[i]);
    if(ans == -1) return {-1, -1};
  }
  return {ans, mod};
}

ll fact[1000010], ifact[1000010];
void makeFac(ll p, ll q) {
  ll pr = 1;
  REP(i, q) pr *= p;
  fact[0] = ifact[0] = 1;
  FOR(i, 1, pr+1) {
    if(i%p == 0) {
      fact[i] = fact[i-1];
    } else {
      fact[i] = fact[i-1] * i % pr;
    }
    ifact[i] = inv(fact[i], pr);
  }
}

ll C(ll n, ll r, ll p, ll q) {
  if(n < 0 || r < 0 || r > n) return 0;
  // pr = p^q
  int pr = 1;
  REP(i, q) pr *= p;

  ll z = n-r;
  int e0 = 0;
  for(ll u=n/p;u;u/=p) e0 += u;
  for(ll u=r/p;u;u/=p) e0 -= u;
  for(ll u=z/p;u;u/=p) e0 -= u;
  int em = 0;
  for(ll u=n/pr;u;u/=p) em += u;
  for(ll u=r/pr;u;u/=p) em -= u;
  for(ll u=z/pr;u;u/=p) em -= u;

  ll ret = 1;
  while(n > 0) {
    ret = ret*fact[n%pr]%pr*ifact[r%pr]%pr*ifact[z%pr]%pr;
    n /= p; r /= p; z /= p;
  }
  (ret *= binpow(p, e0, pr)) %= pr;
  if(!(p==2 && q >= 3) && em%2) ret = (pr-ret) % pr;
  return ret;
}

ll func(ll n, ll r, ll mod) {
  ll x = mod;
  vector<ll> a, m;
  FOR(i, 2, mod+1) if(x%i == 0) {
    ll cnt=0, pr=1;
    while(x%i==0) x/=i, cnt++, pr*=i;
    makeFac(i, cnt);
    a.PB(C(n, r, i, cnt));
    m.PB(pr);
  }

  return crt(a, m).first;
}

signed main(void)
{
  ll n, m, mo;
  cin >> n >> m >> mo;
  cout << func(n, m, mo) << endl;

  return 0;
}

第4回 ドワンゴからの挑戦状 予選 C - Kill/Death

問題ページ C: Kill/Death - 第4回 ドワンゴからの挑戦状 予選 | AtCoder

考えたこと

  • sample1を読むと5デスを4人に割り振るようなのを求められればよさそうでこれは分割数なのでできる
  • sample2を読むとキル数が違う人ごとにまとめてそれぞれ計算していくとよさそう
  • dp[i][j] = (i番目のキル数の人まででjデスを振り分けた場合) みたいなDPが生える
  • 状態数O(Nsum(kill))で遷移O(N)くらいで書けそうなので書くと通った

解法

dp[i][j] = (i番目のキル数の人まででjデスを振り分けた場合) としてDPをする。dp[i][j] = sum(dp[i-1][k] * (j-kデスをi番目のキル数の人数に割り振る)) とするとO(N^2 sum(kill)) で10^7くらいなので通る。j-kデスをi番目のキル数の人数に割り振るパターン数は分割数のDPをしておくと前計算できる。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
#define int ll
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define PB push_back

const ll LLINF = (1LL<<60);
const int INF = (1LL<<30);
const int MOD = 1000000007;

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }
template<class S,class T>
ostream &operator <<(ostream& out,const pair<S,T>& a){
  out<<'('<<a.first<<','<<a.second<<')';
  return out;
}

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

int a[105], b[105];
int dp[1005][1005], dp1[105][1005], dp2[105][1005];
signed main(void)
{
  int n, m;
  cin >> n >> m;
  REP(i, n) cin >> a[i];
  REP(i, m) cin >> b[i];

  VI v1, v2;
  int cnt = 1, sum1 = a[0];
  FOR(i, 1, n) {
    sum1 += a[i];
    if(a[i] == a[i-1]) cnt++;
    else {
      v1.PB(cnt);
      cnt = 1;
    }
  }
  v1.PB(cnt);

  cnt = 1;
  int sum2 = b[0];
  FOR(i, 1, m) {
    sum2 += b[i];
    if(b[i] == b[i-1]) cnt++;
    else {
      v2.PB(cnt);
      cnt = 1;
    }
  }
  v2.PB(cnt);

  // 分割数
  dp[0][0] = 1;
  FOR(i, 1, 1005) {
    FOR(j, 0, 1005) {
      if(j-i >= 0) {
        dp[i][j] = (dp[i-1][j] + dp[i][j-i]) % MOD;
      } else {
        dp[i][j] = dp[i-1][j];
      }
    }
  }

  // v1[0]について
  dp1[0][0] = 1;
  FOR(i, 1, sum2+1) {
    dp1[0][i] = dp[v1[0]][i];
  }

  FOR(i, 1, v1.size()) REP(j, sum2+1) {
    // dp[i-1][k]→dp[i][j]
    REP(k, j+1) {
      // dp[i-1][k] * (j-kをv1[i]人に振り分けたとき)
      (dp1[i][j] += dp1[i-1][k] * dp[v1[i]][j-k] % MOD) %= MOD;
    }
  }

  dp2[0][0] = 1;
  FOR(i, 1, sum1+1) {
    dp2[0][i] = dp[v2[0]][i];
  }
  FOR(i, 1, v2.size()) REP(j, sum1+1) {
    // dp[i-1][k]→dp[i][j]
    REP(k, j+1) {
      // dp[i-1][k] * (j-kをv1[i]人に振り分けたとき)
      (dp2[i][j] += dp2[i-1][k] * dp[v2[i]][j-k] % MOD) %= MOD;
    }
  }

  cout << dp1[v1.size()-1][sum2] * dp2[v2.size()-1][sum1] % MOD << endl;

  return 0;
}

SRM 636 div1 easy ChocolateDividingEasy

解法

二次元累積和を計算しておくと矩形の範囲の値の総和をO(1)で計算できる。あとは切る位置を全通り試すO(H^2W^2)をすればよい。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

int d[55][55] = {};
class ChocolateDividingEasy {
   public:
   int findBest(vector <string> choco)
  {
    int h = choco.size(), w = choco[0].size();
    REP(i, h) REP(j, w) d[i+1][j+1] = choco[i][j]-'0';
    FOR(i, 1, h+1) FOR(j, 1, w+1) d[i][j] += d[i-1][j] + d[i][j-1] - d[i-1][j-1];

    int ret = 0;
    FOR(i, 1, h) FOR(j, i+1, h) FOR(k, 1, w) FOR(l, k+1, w) {
      int v1 = d[i][k],
          v2 = d[i][l] - v1,
          v3 = d[j][k] - v1,
          v4 = d[j][l] - v2 - v3 - v1,
          v5 = d[i][w] - v1 - v2,
          v6 = d[j][w] - v1 - v2 - v3 - v4 - v5,
          v7 = d[h][k] - v1 - v3,
          v8 = d[h][l] - v1 - v2 - v3 - v4 - v7,
          v9 = d[h][w] -v1-v2-v3-v4-v5-v6-v7-v8;
      chmax(ret, min({v1, v2, v3, v4, v5, v6, v7, v8, v9}));
    }

    return ret;
  }
};

SRM635 div1 easy SimilarRatingGraph

解法

i,jから始まる部分列について同型である直線の長さを考える。各区間について傾きが等しく、比率が等しいならば同型であると判定できる。
i+k→i+k+1、j+k→j+k+1の線分について傾きが等しく、day[i+k+1]-day[i+k] と day[j+k+1] - day[j+k] の比率が等しければ部分列[i,i+k]と[j,j+k]は同型である。同型であるところまでkを0~nまで全て試す。
i,jについて全て試す部分でO(N^2)、どこまで同型なのかkを試す部分がO(N)で全体でO(N^3)で求められる。

double = int / int というポカで落として反省。0除算とか色々罠が多そう。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

class SimilarRatingGraph {
   public:
   double maxLength(vector <int> date, vector <int> rating)
  {
    int n = date.size();
    double ret = 0.0;
    REP(i, n) FOR(j, i+1, n) {
      double ratio, len = 0;
      REP(k, n) {
        if(j+k+1 >= n) break;
        // i+k→i+k+1 と j+k→j+k+1 の比が等しいか?
        if((rating[i+k+1]-rating[i+k])*(date[j+k+1]-date[j+k]) == (rating[j+k+1]-rating[j+k])*(date[i+k+1]-date[i+k])) {
          if(k == 0) {
            ratio = (double)(date[j+k+1]-date[j+k])/(date[i+k+1]-date[i+k]);
          } else {
            if(ratio != (double)(date[j+k+1]-date[j+k])/(date[i+k+1]-date[i+k])) {
              break;
            }
          }
          if(ratio > 1) {
            len += sqrt(pow(rating[j+k+1]-rating[j+k], 2) + pow(date[j+k+1]-date[j+k], 2));
          } else {
            len += sqrt(pow(rating[i+k+1]-rating[i+k], 2) + pow(date[i+k+1]-date[i+k], 2));
          }
          chmax(ret, len);
        } else {
          break;
        }
      }
    }

    return ret;
  }
};

SRM634 div1 easy ShoppingSurveyDiv1

概要

M種類の品物を売っている店でN人の人が買い物に来た。同じものを2つ以上買った人はいない。商品iが売れた数はS[i]である。このとき、K個以上の商品を買った人を大口顧客と呼ぶ。このとき、大口顧客の最少人数を求めろ。

解法

できるだけ大口顧客の数を減らしたいのだから、大口顧客は存在する商品を全部買うとした方がよい。あとは大口顧客の数を総当りし、残った商品を購入することが可能かどうか判定し最も少ないものを出力すればよい。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

class ShoppingSurveyDiv1 {
   public:
   int minValue(int N, int K, vector <int> s)
  {
    VI ini = s;
    REP(i, N+1) {
      s = ini;
      int su = 0;
      REP(j, s.size()) {
        s[j] -= i;
        if(s[j] < 0) s[j] = 0;
        su += s[j];
      }
      if((N-i)*(K-1) >= su) {
        return i;
      }
    }

    assert(false);
  }
};

SRM632 div1 easy PotentialArithmeticSequence

概要

正の整数列であるaがある。aは与えられず、d[i] = (a[i]の末尾の0が連続する数) である数列dが与えられる。この数列の部分列で、a[i]が1ずつ増加していく部分列はいくつ存在するか求めろ。

解法

2^6 = 64、n <= 50 から部分列においてd[i] = 6が二つ以上存在することはありえない。したがってd[i]は6以下であるとしてよい。区間[s,f]の部分列について条件をみたすかどうか確認することを考える。d[i]の最大は6なのでmod 64で周期がある。a[s]を1から64として決めおき、区間の要素が全て条件を満たしているかどうか確認すればよい。 計算量は全体で O(2^(max(d[i])) n^3) で解ける。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

class PotentialArithmeticSequence {
   public:
   int numberOfSubsequences(vector <int> d)
  {
    int n = d.size();
    REP(i, n) d[i] = min(d[i], 6);

    // cnt[i] = (iのtrailing zerosの数)
    int cnt[200] = {};
    FOR(i, 1, 200) {
      REP(j, 10) {
        if(i&1<<j) {
          cnt[i] = j;
          break;
        }
      }
    }

    int ret = 0;
    // [i, j] がokかどうか?
    REP(i, n) FOR(j, i, n) {
      // a[i] = sと仮定
      FOR(s, 1, 65) {
        bool flag = true;
        REP(k, j-i+1) {
          // a[i+k] = s+k であるとしたとき問題ないか
          // s+k の trailing zeros の数 == d[i+k] であればok
          if(cnt[s+k] != d[i+k]) {flag = false; break;}
        }
        if(flag) {
          ret++;
          break;
        }
      }
    }

    return ret;
  }
};

SRM631 div1 easy TaroJiroGrid

概要

各マスが黒か白に塗られているn*nの二次元グリッドが与えられる。同じ色がN/2マスより長く連続している列が存在する時、そのグリッドはだめなグリッドである。操作を繰り返すことでだめなグリッドを解消したい。このとき最小の操作回数を求めろ。なお、1回の操作では以下のいずれかを行う。

  • 行を一つ選び黒いマスを全て白に塗る
  • 行を一つ選び白いマスを全て黒に塗る

解法

各操作はある行を全て黒いマスで塗るか白いマスで塗ることと等しい。
N/2マス以上同じ色が存在しないようにするには中央の2行を黒と白でそれぞれ塗ればいいので操作回数は2以下で抑えられる。したがって、操作回数0または1で条件を満たすものがあるか全て試せばよい。愚直に計算してもO(N^4)なので解ける。

#include <bits/stdc++.h>

using namespace std;
typedef long long ll;
typedef vector<int> VI;
typedef vector<VI> VVI;
typedef vector<ll> VL;
typedef vector<VL> VVL;
typedef pair<int, int> PII;

#define FOR(i, a, n) for (ll i = (ll)a; i < (ll)n; ++i)
#define REP(i, n) FOR(i, 0, n)
#define ALL(x) x.begin(), x.end()
#define IN(a, b, x) (a<=x&&x<b)
#define MP make_pair
#define PB push_back
const int INF = (1LL<<30);
const ll LLINF = (1LL<<60);
const double PI = 3.14159265359;
const double EPS = 1e-12;
const int MOD = 1000000007;
//#define int ll

template <typename T> T &chmin(T &a, const T &b) { return a = min(a, b); }
template <typename T> T &chmax(T &a, const T &b) { return a = max(a, b); }

int dx[] = {0, 1, 0, -1}, dy[] = {1, 0, -1, 0};

class TaroJiroGrid {
   public:
   int getNumber(vector <string> g)
  {
    vector<string> s = g;
    int n = g.size();
    int ret = 2;

    REP(i, 2*n+1) {
      g = s;
      int cnt = 0;
      if(i != 0) {
        if(i <= n) g[i-1] = string(n, 'W');
        else g[i-n-1] = string(n, 'B');
        cnt++;
      }

      int ma = 0, con = 1;
      REP(x, n) {
        con = 1;
        FOR(y, 1, n) {
          if(g[y][x] == g[y-1][x]) {
            con++;
          } else {
            ma = max(ma, con);
            con = 1;
          }
        }
        ma = max(ma, con);
      }
      if(ma <= n/2) {
        chmin(ret, cnt);
      }
    }

    return ret;
  }
};