2015-05-23

Codeforces Round #304 (Div. 2) D Soldier and Number Game

Codeforces

dp[i]をi以下の各自然数について素因数分解したときの素因数の個数の総和とすると、dp[a]-dp[b]が答えとなる。よって、前処理で素因数分解をしておけばいい。

エラトステネスの篩をするときにsieve[i]にiの素因数のひとつ(最小や最大のものに固定することもできる)を覚えておくと、素因数分解が簡単にできることを教えてもらった。

#include <bits/stdc++.h>
using namespace std;

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

int sieve[6000000]={};
int dp[6000000]={};

int main(){
    sieve[0]=sieve[1]=1;
    for(int i=2;i*i<=5000000;i++){
        for(int j=2;i*j<=5000000;j++){
            if(sieve[i*j]) continue;
            sieve[i*j]=i;
        }
    }
    dp[0]=dp[1]=0;
	for(int i=2;i<=5000000;i++){
	    if(sieve[i]==0) dp[i]=1;
	    else dp[i] = dp[i/sieve[i]] + 1;
	}
	REP(i,5000000) dp[i+1]+=dp[i];
	int t;
    cin >> t;
    REP(i,t){
        int a,b;
        scanf("%d%d",&a,&b);
	    printf("%d\n",dp[a]-dp[b]);
    }
	return 0;
}

cin,cout は遅すぎ

2015-05-14

高速メビウス変換について

アルゴリズム

高速メビウス変換について自分が勉強した内容をまとめておく。証明はいい感じのサイトを参照のこと。
(http://d.hatena.ne.jp/simezi_tan/20130522/1369203086)

集合に関する関数f,gについて次式が成り立つとき、g(S)が計算できているならば、任意のSについてO(n2^n)でf(S)を計算することができる。
${ \displaystyle f(S) = \sum_{T \subset S} (-1)^{|S| - |T|} g(T) }$

例えば、S={1,2,...,n}として、各i∈Sについて集合A_iが与えられているとき、和集合の要素数を計算したいとする。そして、
${ \displaystyle g(S) = |\bigcap_{i \subset S} A_i| }$
は任意のSについて計算が出来たとする。
このとき、
${ \displaystyle f(S) = \sum_{T \subset S} (-1)^{|S| - |T|} |\bigcap_{i \subset T} A_i| }$
に対して、包除原理の式
${ \displaystyle |\bigcup_{i \in S} A_i| = \sum_{T \subset S} (-1)^{|T|+1} |\bigcap_{j \in T} A_j| }$
を用いると、
${ \displaystyle f(S) = (-1)^{|S| - 1} |\bigcup_{i \subset S} A_i }$
となる。これの絶対値を取れば和集合の要素数を得られる。

さて、f(S)は次の漸化式を考えるとf(S)=f_n(S)として計算できることが証明できる。
${ \displaystyle f_0(S) = g(S) \\ f_k(S) = f_{k-1}(S) , (k \not\in S) \\ f_k(S) = -f_{k-1}(S-\{k\}) + f_{k-1}(S) , (k \in S) }$

bitDPをすることでこの漸化式はO(n2^n)で計算することができる。

これをそのまま書けば次のようになるはず。

//f[n+1][1<<n]
//f[0][T](=g[T]) は計算済みとする
for(int i=1;i<=n;i++)for(int S=0;S<(1<<n);S++){
    if(S & (1<<i)) f[i][S]=-f[i-1][S^(1<<i)]+f[i-1][S];
    else f[i][S]=f[i-1][S];
}

上のDPでは各iにおいて既に計算したSの値しか参照していないため、一次元配列を使いまわすことができる。

for(int i=0;i<n;i++)for(int S=0;S<(1<<n);S++) if(S & (1<<i)) f[S]-=f[S ^ (1<<i)];

2015-04-30

ネットワークフロー

アルゴリズム

練習会でフロー関連の発表をした。内容は↓のとおり。

今までフローを避けてきたので、まとまった勉強をする良い機会となった。しかし、一般マッチングや最大マッチングと最小点カバーにおける関係などはまだ理解できてないので、もう少し蟻本を睨む必要がある。

2015-04-28

SRM 657 Div1 Easy Problem Sets

Topcoder

二分探索するだけ。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
#define ll long long

class ProblemSets {
public:

bool judge(ll E,ll EM,ll M,ll MH,ll H,ll num){
    ll x=max(0ll,num-E);
    ll y=max(0ll,num-H);
    if(x>EM || y>MH) return false;
    if(num>=M && M+EM+MH-x-y<num) return false;;
    return true;
}

long long maxSets(long long E, long long EM, long long M, long long MH, long long H) {
	 ll lb=0;
     ll ub=E+EM+M+MH+H;
     while(ub-lb>1){
        ll mid=(lb+ub)/2;
        if(!judge(E,EM,M,MH,H,mid)) ub=mid;
        else lb=mid;
     }
     return lb;
}
}

2015-04-18

ARC 37 C問題　億マス計算

Atcoder

制約からして二分探索なんだろうとは思ったが本番では解けなかった。
a,bを最初にソートしておき、付値関数C(x)として次のものを考える。

C(x): a[i]*b[j]<=xとなる組(i,j)がK個以上あればtrue、なければfalseを返す

a[i]*b[pos]<=x , a[i]*b[pos+1]>x のとき、a[i]x なるpos'について
pos>pos' が成り立つ。従ってC(x)はO(N)で計算できる。

#include <bits/stdc++.h>
using namespace std;

#define REP(i,n) for(int i=0;i<(int)(n);i++)
typedef long long ll;

bool C(vector<ll> &a,vector<ll> &b,int K,ll x){
    int n=a.size();
    ll res=0;
    int pos=n-1;
    REP(i,n){
        while(pos>=0 && b[pos]>x/a[i]){
            pos--;
        }
        res+=pos+1;
    }
    if(res>=K) return true;
    else return false;
}

int main(){
	int N,K;
    cin >> N >> K;
    vector<ll> a(N),b(N);
    REP(i,N) cin >> a[i];
    REP(i,N) cin >> b[i];
    sort(a.begin(),a.end());
    sort(b.begin(),b.end());
    ll lb=0,ub=a[N-1]*b[N-1]+2;
    while(ub-lb>1){
        ll mid=(lb+ub)/2;
        if(C(a,b,K,mid)) ub=mid;
        else lb=mid;
    }
    cout << ub << endl;
	return 0;
}

2015-04-18

最短経路問題

アルゴリズム

蟻本等の内容を復習していく。

Bellman-Ford法

頂点sから頂点iまでの仮の最短距離をd[i]とする
0: d[s]=0, d[i(≠s)]=INF とする
1: u,v∈V に対してd[v]=min(d[v],d[u]+w(u,v))で更新。更新できなくなったら終了
step1がV回実行されても終了しない場合、負の閉路が存在することが分かる。

計算量はstep1の操作がO(V)で高々E回実行されるからO(VE)

//Bellman-Ford algorithm
struct edge{int from,to,cost;};
const int INF=1e9;

bool BellmanFord(vector<vector<edge> > &g,vector<int> &d,vector<int> &prev,int s){
    int V=g.size();
    for(int i=0;i<V;i++) d[i]=(i==s)?0:INF;
    for(int i=0;i<V;i++){
        bool update=false;
        for(int u=0;u<V;u++){
            for(int j=0;j<g[u].size();j++){
                edge e=g[u][j];
                if(d[e.to]>d[e.from]+e.cost){
                    update=true;
                    d[e.to]=d[e.from]+e.cost;
                    prev[e.to]=e.from;
                }
            }
        }
        if(!update) return true;
    }
    return false;   //there is a negative loop
}
//
vector<int> buildPath(vector<int> &prev,int t){
    vector<int> res;
    for(int v=t;v>=0;v=prev[v]){
        res.push_back(v);
    }
    reverse(res.begin(),res.end());
    return res;
}

Dijkstra法

頂点sから頂点iまでの仮の最短距離をd[i]とする。Bellman-Ford法では次に使う頂点の探索に無駄があったが、Dijlkstra法では(仮の最短距離、頂点)の候補をヒープで管理することで高速で計算できる。但し、負の辺が含まれる場合は使うことができない。

計算量は、次に使う頂点の探索がO(logV)で全体ではO(E logV)

//Dijkstra
struct edge{int from,to,cost;};
const int INF=1e9;

void Dijkstra(vector<vector<edge> > &g,vector<int> &d,int s){
    int V=g.size();
    for(int i=0;i<V;i++) d[i]=(i==s)?0:INF;
    priority_queue<pair<int,int>,vector<pair<int,int> >,greater<pair<int,int> > > que;
    que.push(make_pair(0,s));
    while(!que.empty()){
        pair<int,int> p=que.top();que.pop();
        int u=p.second;
        if(d[u]<p.first) continue;
        for(int i=0;i<g[u].size();i++){
            edge e=g[u][i];
            if(d[e.to]>d[e.from]+e.cost){
                d[e.to]=d[e.from]+e.cost;
                que.push(make_pair(d[e.to],e.to));
            }
        }
    }
}

2015-04-17

SRM 656 Div1 Easy RandomPancakeStack

Topcoder

dp(i,j):i番目の操作でj番目のパンケーキを置く確率

class RandomPancakeStack {
public:

int N;
double memo[300][300];

double dp(int i,int j){
    if(memo[i][j]>1e-9) return memo[i][j];
    if(i==0) return memo[i][j]=1.0/N;
    double res=0;
    for(int k=j+1;k<N;k++) res+=dp(i-1,k)/(N-i);
    return memo[i][j]=res;
}

double expectedDeliciousness(vector <int> d) {
	N=d.size();
    REP(i,300)REP(j,300) memo[i][j]=0;
    double res=0;
    REP(i,N)REP(j,N) res+=dp(i,j)*d[j];
    return res;
}