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;
}

2015-04-17

AOJ 1056 Ben Toh

AOJ

n日までに得る弁当の総数の期待値を計算すればいいわけだが、漸化式がうまく作れない && 他の人のコードを読んでも自明のようにdpしていてよく分からない、というツラさがあったので、この問題の漸化式が私のように漸化式が思いつかなかった人のために愚直に式の導出をする(もっとスマートに立式できる方法があったら教えて下さい)。

$E_n$ : n日までに得る弁当の総数
$F_n$ : n日目に弁当を得る確率
$G_{n,k}$ : n日目に弁当を得て、それまでにk-1回連続で弁当を得ている確率

期待値の性質から $E_n=\sum_{i=1}^n F_i$ なので、 $F_n$ の漸化式を求めればよい。

ここで $G_{n,k}(1 \leq k \leq n)$ について以下の漸化式が成り立つ。
$\displaystyle G_{n,k}=G_{n-1,k-1} \times \frac{1}{2^{k-1}} \ (1 \leq k < n) \\ G_{n,1}=1-F_n$
上式を繰り返し用いて、
$\displaystyle G_{n,k}=...=G_{n-k+1,1} \times \prod_{i=1}^{k-1} \frac{1}{2^i}=G_{n-k+1,1} \times 2^{-\frac{k(k-1)}{2}}$
ゆえに、
$\displaystyle G_{n,k} = \begin{cases} (1-F_{n-k})2^{-\frac{k(k-1)}{2}} & (1 \leq k < n) \\ 2^{-\frac{n(n-1)}{2}} & (k=n) \end{cases}$
ここで、 $F_{n}=\sum_{i=1}^n G_{n,i}$ であるから $n \geq 2$ のとき、
$\displaystyle F_n = 2^{-\frac{n(n-1)}{2}} + \sum_{i=1}^{n-1}{(1-F_{n-k}) 2^{-\frac{k(k-1)}{2}}}$
$F_1=1$ であるから、これで $F_n$ を逐次的に求めることができる。

解答の上では $2^{\frac{k(k-1)}{2}}$ がdouble型の有効桁数よりも小さい範囲で和を取れば十分である。

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

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

double memo[100001]={};

double cff(int k){
    if(k>=20) return 0.0;
    return pow(2.0,-k*(k-1)/2);
}

double dp(int n){
    if(n==1) return memo[1]=1.0;
    if(memo[n]>1e-10) return memo[n];
    double res=cff(n);
    for(int k=1;k<=min(n-1,20);k++){
        res+=(1-dp(n-k))*cff(k);
    }
    return memo[n]=res;
}

int main(){
	int n;
    while(cin >> n && n){
        double res=0;
        for(int i=1;i<=n;i++) res+=dp(i);
        cout << fixed << setprecision(10) << res << endl;
    }
	return 0;
}

2015-04-16

scipyの特殊関数

python

量子力学の教科書を読んでいて、球Bessel関数とかLaguerreの陪多項式がある微分方程式の解だと言われてもいまいちピンと来なかったので適当に実装してmatplotlibでグラフを書いてみようと思っていたところ、scipyに特殊関数のmoduleがあることを知った。

Special functions (scipy.special) — SciPy v0.15.1 Reference Guide

有名所の特殊関数はだいたい実装されており、折角なので特殊関数の概形を見てみるだけでなく水素原子のSchrödinger方程式の解もプロットしてみることにした。

ここで問題なのは3次元上の各点での波動関数の振幅の2乗をカラースケールで表示しても、意味の分からないグラフとなってしまうことだ。これを回避するために、各座標に点がプロットされる確率を波動関数の振幅の２乗と対応するようにした。n=1,2ではこれでうまくいくが、n>=3だと空間的に広がるために毎回パラメターをいじる必要がある。

下図は(n,l,m)=(2,1,0)のときの結果。ボーア半径を1としてある。
f:id:lan496:20150416001412p:plain

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import scipy.special
import numpy as np
import random
import math

def freeParticle(n,l,m,r,theta,phi):
    res = (-1)**(m+abs(m))
    res *= scipy.special.sph_harm(m,l,phi,theta)
    res *= scipy.special.sph_jn(0,r)[0][0]
    return res

def hydrogen(n,l,m,r,theta,phi):
    res = (-1)**(m+abs(m))
    res *= scipy.special.sph_harm(m,l,phi,theta)
    
    rho = 2.0*r/n
    res *= scipy.special.assoc_laguerre(rho,n+l,2*l+1)
    tmp = np.sqrt(math.factorial(n-l-1)/2.0/n*math.factorial(n+l)**-3*(2.0/n)**3)
    res *= -tmp
    res *= np.exp(-rho/2.0)*(rho**l)
    
    return res

def plotGraph(n,l,m,waveFunction):
    L = 2.0
    dl = 0.1

    fig = plt.figure()
    ax = Axes3D(fig)
    
    ax.set_xlabel("X-axis")
    ax.set_ylabel("Y-axis")
    ax.set_zlabel("Z-axis")
    
    ax.set_xlim(-L,L)
    ax.set_ylim(-L,L)
    ax.set_zlim(-L,L)

    for z in np.arange(-L,L,dl):
        for y in np.arange(-L,L,dl):
            for x in np.arange(-L,L,dl):
                r = np.sqrt(x**2 + y**2 + z**2)
                theta = np.arccos(z/r)
                phi = np.arctan2(y,x)

                func = waveFunction(n,l,m,r,theta,phi)
                probability = func*np.conj(func)
                print probability
                ran = [random.random() for i in range(20)]
                if min(ran) < probability:
                    ax.scatter3D(x,y,z,s=10)

    plt.show()

if __name__=='__main__':
    plotGraph(2,1,-1,hydrogen)

2015-04-15

AOJ 1163 Cards

AOJ

Ford-Fulkerson法を最近勉強したのでやっと解けるようになった。

二部グラフの最大マッチングのサイズを求めるという典型題中の典型題。

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

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

//Ford-Fulkerson's algorithm
struct edge{int to,cap,rev;};
const int INF=1e9;

void addEdge(vector<vector<edge> > &g,int from,int to,int cap){
    g[from].push_back((edge){to,cap,(int)g[to].size()});
    g[to].push_back((edge){from,0,(int)g[from].size()-1});
}

int dfs(vector<vector<edge> > &g,vector<bool> &used,int v,int t,int f){
    if(v==t) return f;
    used[v]=true;
    for(int i=0;i<(int)g[v].size();i++){
        edge& e=g[v][i];
        if(!used[e.to] && e.cap>0){
            int d=dfs(g,used,e.to,t,min(f,e.cap));
            if(d>0){
                e.cap-=d;
                g[e.to][e.rev].cap+=d;
                return d;
            }
        }
    }
    return 0;
}

int FordFulkerson(vector<vector<edge> > &g,int s,int t){
    int flow=0;
    for(;;){
        vector<bool> used(g.size(),false);
        int f=dfs(g,used,s,t,INF);
        if(f==0) return flow;
        flow+=f;
    }
}
//

int gcd(int a,int b){
    if(a<b) swap(a,b);
    if(b==0) return a;
    return gcd(b,a%b);
}

int main(){
	int m,n;
    while(cin >> m >> n && m){
        vector<int> b(m),r(n);
        REP(i,m) cin >> b[i];
        REP(i,n) cin >> r[i];
        vector<vector<edge> > g(m+n+2);
        REP(i,m) addEdge(g,0,i+1,1);
        REP(i,n) addEdge(g,m+1+i,m+n+1,1);
        REP(i,m)REP(j,n){
            if(gcd(b[i],r[j])>1) addEdge(g,i+1,m+1+j,1);
        }
        cout << FordFulkerson(g,0,m+n+1) << endl;
    }
	return 0;
}

2015-04-12

ABC 31 D 多重ループ

Atcoder

$\begin{eqnarray}{}_n H _k = _{n-1+k} C _{k-1} \end{eqnarray}$ をmod 1e9+7 での逆元を使いながら計算すれば良い。

最初 #define MOD 1000000007 としていたが、これだとMODはint型と判断されオーバーフローを起こしてしまってつらかった。

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

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

const ll MOD=1000000007;

ll inverse(ll a){
    ll res=1;
    ll base=a;
    REP(i,33){
        if(((MOD-2)>>i) & 1) res=(res*base)%MOD;
        base=(base*base)%MOD;
    }
    return res;
}

ll fact(ll a){
    ll res=1;
    for(int i=1;i<=a;i++) res=(res*i)%MOD;
    return res;
}

ll nHk(ll n, ll k){
    ll a=fact(n-1+k);
    ll b=inverse(fact(k));
    ll c=inverse(fact(n-1));
    return (((a*b)%MOD)*c)%MOD;
}

int main(){
	int n,k;
    cin >> n >> k;
    cout << nHk(n,k) << endl;
	return 0;
}

2015-04-11

AOJ 2587 Broken Cipher Generator

AOJ

以下のBNFに従ってparserを書く。

・<Cipher>:=(<String>)+
・<String>:='[' + <Cipher> + ']' | <Letter>
・<Letter>:=('+'|'-')* + <Alphabet>

但し、'?'を含むの要素は'A'と同一視する。

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

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

string expr_string();
char expr_letter();

int c;
string s;

string expr_cipher(){
    string res="";
    while(c!=s.length() && s[c]!=']'){
        res+=expr_string();
    }
    return res;
}

string expr_string(){
    string res="";
    char p=s[c];
    if(p=='['){
        c++;
        res=expr_cipher();
        reverse(res.begin(),res.end());
        c++;
    }else{
        res+=expr_letter();
    }
    return res;
}

char expr_letter(){
    int cnt=0;
    while(s[c]=='+' || s[c]=='-'){
        if(s[c]=='+') cnt++;
        else cnt--;
        c++;
    }
    if(s[c]=='?'){
        c++;
        return 'A';
    }
    char res='A'+(s[c]-'A'+cnt+130)%26;
    c++;
    return res;
}

int main(){
    while(cin >> s && s!="."){
        c=0;
        cout << expr_cipher() << endl;
    }
	return 0;
}