Richardson補外とか丸め誤差について

数値計算の常識』(伊理正夫・藤野和建著)を読んでいろいろ試したくなったので、取り敢えず簡単にできるものをやってみる。

積分\( I = \int_{a}^{b} f(x) dx\)をN等分した台形\(I_N\)で近似するとき、
\( I_{1} = \frac{1}{2}(b - a)(f(a) + f(b)) \)
\( I_{2^{N+1}} = \frac{1}{2} (I_{2^N} + h \sum_{i = 1}^{2^N} f(a + (i - \frac{1}{2})h) ) \), 但し\( h = \frac{b - a}{2^N}\)
で近似値を(逐次的に)出すことができる。

また、f(x)が\(C^{2p}\)級のとき、hに関係しない係数\( d_{2j}\)を用いて
(*)\( I_N - I = \sum_{i=1}^{p} d_{2i} h^{2i} + O(h^{2p+2})\)
と表すことができるらしい(間違っていたらすいません)。

(*)より、上手く数値積分ができているならば、k=1,...,pとして
\( I_{N}^{k} = \frac{ I_{N}^{k - 1} - 4^{-k} I_{N / 2}^{k - 1}}{ 1 - 4^{-k}} \)
但し、\( I_{N}^{0} = I_{N}\)と定めたとき[Richardson補外]、少し計算すると\( I_{N}^{k} - N\)はNが2倍になると\(1 / 4^{k+1}\)倍になること、即ち\( I_{N}^{k} - I_{N/2}^{k}\)はNが2倍になると\(1 / 4^{k+1}\)倍になることが分かる。これより、もしも\( I_{N}^{k}\)の階差を対数目盛でみたときに直線にのっていなければ丸め誤差等の影響で健全な計算が行えていないことが分かる。ちなみに、k=1の場合に得られる近似式はSimpsonの式に一致する。


例として次の(よく例題に出てくる)積分数値計算で求めることを考える。
\( I = \int_{-1}^{1} \frac{2}{1 + x^2} dx = \pi \)

2^34(=17179869184)等分まで倍精度で求めた。更に、ここから\( I_{N}^{1} \)と \(I_{N}^{2}\)を計算した。

\(I_{N}^{k}\)階差をプロットしたものが下図(順に\(I_{N}, I_{N}^{1}, I_{N}^{2}\)の場合)。

f:id:lan496:20150906011606p:plain
\( N \le 2^{20} \) くらいまでは順調そうだがそれより大きくなると丸め誤差が無視できなくなるようだ。\(I_{2^{20}} = 3.141592653589158068\)であった。

f:id:lan496:20150906011615p:plain
\( N \le 2^9 \) くらいでもう怪しくなっている。\( I_{2^9}^{1} = 3.141592653589792228 \)であった。

f:id:lan496:20150906011620p:plain
ダメダメ。


結果、\(I_{N}^{1}\)を活用することで小数点桁数13桁くらい(3.1415926535897...)まで求めることができた。

実際には直線に乗らない範囲で円周率により近い値を得ているが、本来真値を知っている訳ではないのでそれをより真値に近いものと判断する根拠は乏しい。

[参考]

#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <cstdio>
using namespace std;

double a = -1.0;
double b = 1.0;

double func(double x){
  return 2.0 / (1.0 + x * x);
}

int main() {
  int M = 34;
  vector<double> integral(M);
  integral[0] = (b - a) * (func(a) + func(b)) / 2.0;
  
  for(int i=1;i<M;++i){
    double middle = 0.0;
    long long N = pow(2, i - 1);
    double h = (b - a) / N;
    
    for(long long i=1;i<=N;i++){
      middle += h * func(a + (i - 0.5) * h);
    }
    
    integral[i] = (integral[i-1] + middle) / 2.0;
  }
  for(int i=0;i<M;i++){
    printf("%lld %.18lf\n", (long long)pow(2, i), integral[i]);
  }
  return 0;
}
N \(I_{N}\) \(I_{N}^{1}\) \( I_{N}^{2}\)
1 2.000000000000000000 0.000000000000000000 0.000000000000000000
2 3.000000000000000000 3.333333333333333037 0.000000000000000000
4 3.100000000000000089 3.133333333333333304 3.555555555555555358
8 3.131176470588235450 3.141568627450980422 3.119999999999999662
16 3.138988494491089298 3.141592502458707248 3.142117647058823682
32 3.140941612041388886 3.141592651224821786 3.141594094125888859
64 3.141429893174974453 3.141592653552836456 3.141592661142562637
128 3.141551963485655019 3.141592653589214912 3.141592653708037641
256 3.141582481063752041 3.141592653589784234 3.141592653591640083
512 3.141590110458282403 3.141592653589792228 3.141592653589822426
1024 3.141592017806918768 3.141592653589797557 3.141592653589792672
2048 3.141592494644073863 3.141592653589791784 3.141592653589798001
4096 3.141592613853361193 3.141592653589790007 3.141592653589791340
8192 3.141592643655684025 3.141592653589791784 3.141592653589790007
16384 3.141592651106267731 3.141592653589795336 3.141592653589791784
32768 3.141592652968914656 3.141592653589796669 3.141592653589795781
65536 3.141592653434578608 3.141592653589800221 3.141592653589796669
131072 3.141592653550970837 3.141592653589768247 3.141592653589800221
262144 3.141592653580073335 3.141592653589774020 3.141592653589766471
524288 3.141592653587365724 3.141592653589796669 3.141592653589774464
1048576 3.141592653589158068 3.141592653589754924 3.141592653589798001
2097152 3.141592653589570183 3.141592653589707851 3.141592653589752260
4194304 3.141592653589582618 3.141592653589587059 3.141592653589704742
8388608 3.141592653589698081 3.141592653589736273 3.141592653589578621
16777216 3.141592653589746931 3.141592653589763362 3.141592653589746487
33554432 3.141592653589535988 3.141592653589465378 3.141592653589765582
67108864 3.141592653589374340 3.141592653589320605 3.141592653589445838
134217728 3.141592653588935136 3.141592653588788586 3.141592653589311279
268435456 3.141592653589452944 3.141592653589625250 3.141592653588753059
536870912 3.141592653589918793 3.141592653590073780 3.141592653589681206
1073741824 3.141592653590653761 3.141592653590898454 3.141592653590103534
2147483648 3.141592653591862572 3.141592653592265805 3.141592653590953521
4294967296 3.141592653595981055 3.141592653597353735 3.141592653592356843
8589934592 3.141592653606938512 3.141592653610590702 3.141592653597693019

AOJ 2333 My friends are small

まず、wを昇順でソートする。
dp[i][j] = #{w[i]以降のみを使って和がjとなる組}
sum[i] = w[i]までの和
とするとき、次の3つの場合を考えれば良い。

(i) w[0]を使わない
(ii) w[0]からw[i]までを使い、w[i+1]を使わない(0 <= i < N-1)
(iii) 全て使う

(i)の場合, (i=1,j) s.t. 0 <= W - j < w[0] について、dp[i][j]をresに加える
(ii)の場合、 (i,j) s.t. 0 <= W - sum[i] - j < w[i+1] について、dp[i][j]をresに加える
(iii)の場合、 W >= sum[N-1] ならばresに1を加える。

#include <bits/stdc++.h>

#define FOR(i,a,b) for(int i=(a);i<(int)(b);++i)
#define REP(i,n) FOR(i,0,n)
#define ALL(x) (x).begin(),(x).end()

using namespace std;

const int INF = 1e9;
const int MOD = 1000000007;

int main() {
  int N,W;
  cin >> N >> W;
  vector<int> w(N),sum(N);
  REP(i,N) cin >> w[i];
  sort(ALL(w));
  sum[0] = w[0];
  REP(i,N-1) sum[i+1] = sum[i] + w[i+1];
  
  vector<vector<int>> dp(N+1, vector<int> (W+1,0));
  dp[N][0] = 1;
  for(int i=N-1;i>=0;i--)REP(j,W+1){
    if(j >= w[i]) dp[i][j] += dp[i+1][j-w[i]];
    dp[i][j] += dp[i+1][j];
    dp[i][j] %= MOD;
  }

  int res = 0;
  if(W >= sum[N-1]) res++;
  REP(j,W+1) if(j > W - w[0]){
    res = (res + dp[1][j]) % MOD;
  }
  REP(i,N-1)REP(j,W+1){
    if(j > W - sum[i+1] && j <= W - sum[i]){
      res = (res + dp[i+2][j]) % MOD;
    }
  }
  cout << res << endl;
  return 0;
}

SRM 664 Div1 Easy Bear Plays

組(X,Y)(X <=Y)を一回操作して(2X,Y-X)にするとき、石の総量S = X + Y は変化しない。よって、少ない方の個数がXのとき、一回の操作後に得られる組の少ない方の個数f(X)は

f(X) = 2X or S - 2X

これをmod S で見ると、f(X) ≡ 2X or -2X であるから、f^K (X) ≡ 2^K * X or -2^K * X

f^K (X) ≡ -2^K * X は個数がS - 2^K * X の場合に対応するから、X = min(A,B)として、min(2^K * X, S - 2^K * X) (mod S) が答えとなる。

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

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

typedef long long ll;

class BearPlays{
public:
  ll modpow(ll K, ll S){
    if(K == 0)return 1;
    if(K % 2) return (2 * modpow(K - 1, S)) % S;
    return (modpow(K/2, S) * modpow(K/2, S)) % S;
  };
  int pileSize(int A, int B, int K){
    ll S = A + B;
    A = min(A, B);
    ll res = (A * modpow(K,S)) % S;
    return min(res, S - res);
  };
};

Codeforces Round #304 (Div. 2) E Soldier and Traveling

sourceとsinkを追加し流量制限にa[i],b[i]を与えることで二部グラフのフローになることは本番中に気づいたが、YESだった場合に経路復元する方法が分からなかった。

考えてみると、逆辺のcapはその辺に流したフローの量であるから、ここから容易に復元することができる。

それと、source=2n, sink=2n+1 みたいに置いておくと、デバッグしやすいことを教えてもらった。

#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 main(){
	int n,m;
    cin >> n >> m;
    vector<int> a(n),b(n);
    REP(i,n) cin >> a[i];
    REP(i,n) cin >> b[i];
    vector<vector<edge> > g(2*n+2);
    int source=2*n;
    int sink=2*n+1;
    REP(i,n) addEdge(g,source,i,a[i]);
    REP(i,n) addEdge(g,n+i,sink,b[i]);
    REP(i,n) addEdge(g,i,n+i,INF);
    REP(i,m){
        int p,q;
        cin >> p >> q;
        addEdge(g,p-1,q-1+n,INF);
        addEdge(g,q-1,p-1+n,INF);
    }
    int sum1=0,sum2=0;
    REP(i,n){
        sum1+=a[i];
        sum2+=b[i];
    }
    int f=FordFulkerson(g,source,sink);
    if(f!=max(sum1,sum2)){
        cout << "NO" << endl;
    }else{
        vector<vector<int> > res(n,vector<int> (n));
        REP(i,n){
            for(auto e : g[i+n]){
                if(e.to<n) res[e.to][i]=e.cap;
            }
        }
        cout << "YES" << endl;
        REP(i,n){
            REP(j,n){
                if(j) cout << " ";
                cout << res[i][j];
            }
            cout << endl;
        }
    }
	return 0;
}

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

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 は遅すぎ

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

高速メビウス変換について自分が勉強した内容をまとめておく。証明はいい感じのサイトを参照のこと。
(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)];

ネットワークフロー

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

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