Approximation Algorithms Part I Week 2

　Week2 の講義を視聴した。Week2のテーマはナップサック問題だった。

　ナップサック問題は、動的計画法の基礎中の基礎ということで、プロコンで何度も解いているので最初の方はだいたい知っている話だった。価値が小さな整数の場合は多項式時間の動的計画法で解けるが、任意の実数の場合はどうすればよいか？ Scaling、Roundingというテクニックを使うことで近似解を求めることができる。

　Roundingのテクニックは Week1 の vertex cover にも出来てたが、vertex cover の場合は LP relaxation を解いた出力を rounding するのに対して、ナップサック問題の場合は 入力値を rounding する。アイデア自体はシンプルで実数を整数に変換して動的計画で解けばOK。

　Scalingのパラメータの取り方によって、任意のε > 0 に対して|OPT - ApproximationOPT| < εとなるような近似問題を定式化できるが、実行時間は1/εに対して線形で増えていくので実用的には程よいεを見極める必要がある。このあたりの reasoning/analysis は面白かった。

　

Approximation Algorithms Part I Week 1

　Courseraで「Approximation Algorithms Part I」という講義を取り始めた。以下のようなことを教えてくれる。

• クラスNPの組み合わせ最適化問題の近似解を多項式アルゴリズムで求める
• 近似アルゴリズムの精度
• 近似解が実行可能であることの考察

ちなみにコースを提供しているENS（École normale supérieure; 高等師範学校）というのはフランスの大学らしい。

　第1週は Vertex Cover を解くための Linear Programming RelaxationとRoundingについてだった。さっき課題を提出したが、記述式の問題で、生徒同士で評価しあう形式だった。数式を書かないといけないのでひさびさにlatexを使って学生に戻った気分になった。

　組み合わせ最適化問題は実世界の問題にも適用できる場面が多いので、基本を抑えてエンジニアとしての引き出しを増やしておきたい。

関数の括弧の中のセミコロン

　関数の定義で f(x; θ) のような表記があった。

xとθの関数なので f(x, θ) でいいのではと思ったが、セミコロンで区切られている。
セミコロンで区切ると、セミコロンの前は変数、セミコロンの後はパラメータという意味らしい。

f(x; θ)はパラメータθによって特徴付けられる、xを変数として受け取る関数fみたいなイメージ。
θは変数ではなく、fを定義するときに決定する固定値というふうに考える。

集合族とは

　グラフィカルモデルの勉強をしていたら、集合族という聞きなれない言葉に遭遇したので調べてみた。

集合族とは、集合のあつまりのこと。
特に、ある集合Sの部分集合族とは、Sの部分集合のあつまりのことである。

具体的な例を考えてみる。
集合SをS = {1,2,3,4}と定義する。
以下にSの部分集合族をいくつか列挙してみる。

Sのべき集合は、Sの部分集合族である。
べき集合とは、以下のようにSの考えられる部分集合をすべて列挙したものである。
{{}, {1}, {2}, {3}, {4}, {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, {3,4}, {1,2,3}, {1,3,4}, {1,2,4}, {2,3,4}, {1,2,3,4}}

Sのk-元部分集合からなる集合は、Sの部分集合族である。
k-元部分集合とは要素数がkの部分集合のこと。
例えば以下のようなSの1-元部分集合は、Sの部分集合族である。
{{1}, {2}, {3}, {4}}

適当にSから部分集合を選んでそれを要素とする集合を作ると、それはSの部分集合族である。
例えば、{{1,2}, {1,2,3}, {3,4}, {4}}など。

直線型グラフのマッチング数

　以下のような直線型グラフのマッチングについて考える。

一般にノード数nの直線型グラフからk個のマッチングを選ぶパターン数は_n-kC_kとなる。

例えば上のグラフの場合（n=6）は、
k = 1なら_6-1C₁ = 5、
k = 2なら_6-2C₂ = 6、
k = 3なら_6-3C₃ = 1
となり、実際に個数を数えてみると正しいことが分かる。

何をやってるかというと、マッチングに使われる枝の左側のノードを選んでいるわけである。普通にノードを選ぶだけだと一つのノードが複数のマッチングに属してしまうかもしれないので、あらかじめ右側の端点に使われるノードを除いておく。

すると左側のノードの選び方により一意にマッチングが決定する。逆にマッチングの選び方に対して、左側のノードの選び方は一意に決まる。よって一対一の対応にあるのでこの方法で数え上げができることが分かる。

モンモール数を包除原理で計算する

　モンモール数を計算するときは漸化式を使えば計算できるが、実は包除原理を使っても計算できるのではないかと思ったので試してみた。

　予想どおり包除原理でも計算できるようだ。
計算量は漸化式を使うとO(N)、包除原理を使うとO(N²)なので漸化式の方が高速。

#include <iostream>
#include <cassert>

using namespace std;

const long long MOD = 1e9 + 7;
const int N = 1000;
long long f[N];
long long comb[N][N];
long long x[N];
long long y[N];

// 漸化式で解く
void solveWithRec() {
  x[0] = 1, x[1] = 0;
  for (int i = 2; i < N; i++)
    x[i] = (i-1) * (x[i-1] + x[i-2]) % MOD;
}

// 包除原理で解く
void solveWithIncExc() {
  f[0] = 1;
  for (int i = 1; i < N; i++)
    f[i] = i * f[i-1] % MOD;

  comb[0][0] = 1;
  for (int i = 1; i < N; i++) {
    comb[i][0] = comb[i][i] = 1;
    for (int j = 1; j < i; j++)
      comb[i][j] = (comb[i-1][j] + comb[i-1][j-1]) % MOD;
  }
  for (int i = 0; i < N; i++) {
    y[i] = f[i];
    for (int j = 1; j <= i; j++) {
      if (j % 2)
        y[i] -= comb[i][j] * f[i-j] % MOD;
      else
        y[i] += comb[i][j] * f[i-j] % MOD;
      y[i] = (y[i] % MOD + MOD) % MOD;
    }
  }
}

int main(int argc, char *argv[])
{
  solveWithRec();
  solveWithIncExc();

  for (int i = 0; i < N; i++)
    assert(x[i] == y[i]);

  cout << "same!" << endl;
  
  return 0;
}

SRM 507 Div1 500 CubePacking

問題
辺の長さが1の立方体がNs個、辺の長さがLの立方体がNb個ある。
これらの立方体を出来るだけ体積が小さい直方体に格納したい。
直方体の体積の最小値を求めよ。

解法
解の下限はNb*L*L*L+Nsである。
解の上限はNb*L*L*L+Ns+L*L-1である。（L*Lの底面を持つ直方体に立方体を詰め込む）
ということで、 解の候補をすべて列挙しても高々L*L程度のループをまわせばよい。

あとは、それぞれの解の候補に対して、その体積になるような直方体を全列挙すればよい。これは時間がかかりそうだが約数の個数は思ったほど多くないため十分に間に合う。

以下に整数の約数の個数がどの程度かを示す。
この結果はN周辺の数をサンプリングし、それらの約数の個数の平均値を取ったものである。

N	約数の個数
106	14.9
107	17.3
108	19.6
109	21.9

ざっくりと高々log(N)程度と覚えておくと良さそう（あくまでも平均値であることに注意）。

ソースコード

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

vector<long long> divisor(long long n) {
  vector<long long> ret;
  for (long long i = 1; i * i <= n; i++) {
    if (n % i == 0) {
      ret.push_back(i);
      if (i != n / i)
        ret.push_back(n / i);
    }
  }
  return ret;
}

class CubePacking {
  public:
  int getMinimumVolume(int Ns, int Nb, int L) {
    int s = L*L*L*Nb+Ns;
    for (int x = s;; x++) {
      auto ds = divisor(x);
      REP (i, ds.size()) REP (j, i+1) {
        int a = ds[i];
        int b = ds[j];
        int c = x / a / b;
        if (a * b * c != x)
          continue;
        if ((a/L) * (b/L) * (c/L) >= Nb)
          return x;
      }
    }
    return -1;
  }
};

SRM 657 Div2 1000 PolynomialRemainder

問題

a, b, cが与えられる。

a x^2 + b x + c = 0 (mod 1,000,000,000)

となるようなxを求めよ。

解法

10^9で割って0になるので、

2^9で割って0、5^9で割って0となるようなxを求めて、中国剰余定理をすればよい。

実装

using namespace std;

#define ALL(x) (x).begin(), (x).end()
#define EACH(itr,c) for(__typeof((c).begin()) itr=(c).begin(); itr!=(c).end(); itr++)  
#define FOR(i,b,e) for (int i=(int)(b); i<(int)(e); i++)
#define MP(x,y) make_pair(x,y)
#define REP(i,n) for(int i=0; i<(int)(n); i++)

long long extgcd(long long a, long long b, long long &x, long long &y) {
  long long d = a;
  if (b != 0) {
    d = extgcd(b, a%b, y, x);
    y -= a/b * x;
  } else {
    x = 1;
    y = 0;
  }
  return d;
}

pair<long long, long long> crm(long long a, long long n, long long b, long long m) {
  a %= n;
  b %= m;
  long long p, q;
  extgcd(n, m, p, q);
  long long mod = n * m;
  long long x = (a * q  * m) % mod + (b * p * n) % mod;
  x = (x % mod + mod) % mod;
  return make_pair(x, mod);
}

long long calc(long long a, long long b, long long c, long long x, long long mod) {
  long long ret = ((a % mod * x) % mod) * x % mod;
  ret = (ret + b * x) % mod;
  ret = (ret + c) % mod;
  return ret;
}

class PolynomialRemainder {
public:
  int findRoot(int a, int b, int c) {
    long long s, t;

    int n = 1;
    REP (i, 9) n *= 2;
    for (s = 0; s < n; s++) {
      if (calc(a,b,c,s,n) == 0)
        break;
    }
    if (s == n)
      return -1;

    int m = 1;
    REP (i, 9) m *= 5;
    for (t = 0; t < m; t++) {
      if (calc(a,b,c,t,m) == 0)
        break;
    }
    if (t == m)
      return -1;

    cout << s << " " << t << endl;
    auto ret = crm(s, n, t, m);
    return ret.first;
  }
};

Codeforces Round #338 Multipliers

問題概要
とても大きな数nがある。
nの素因数が与えられるので、nの約数の総積をmod 1000000007で求めよ。

制約
nの素因数の個数 <= 200,000
2 <= nの素因数 <= 200,000

解き方
数論の知識が必要。

1. d(x)をxの約数の個数とすると、gcd(a, b) = 1のとき、d(ab) = d(a) d(b)となる。
2. nの約数の総積は、n^(d(n)/2)となる。
3. gcd(a, n) = 1のとき、a^x = a^(x % (n-1)) (mod n)となる。

1. はdivisor関数はmultiplicativeという有名な話。因数分解された数の約数の個数を高速に数える方法を理解していれば確かにそうだと分かる。

2. は最近知った。これは、約数を昇順に列挙したものと、降順に列挙したものを上下に並べて、上下の積を取ればそうなることが分かる。

3. は知らなかった。modで計算するとき冪乗の右肩の数は%取っていけないのは知っていたけど、mod nのときは% (n-1)を計算すればいいらしい。これはフェルマーの小定理から導ける。

サンプル解
上の3つを組み合わせて解く。詳しい式変形はeditorial参照。

#include <iostream>
#include <map>

using namespace std;

const long long MOD = 1e9 + 7;

long long mpow(long long x, long long p, long long mod) {
  long long ret = 1;
  while (p) {
    if (p & 1)
      ret = ret * x % mod;
    x = x * x % mod;
    p >>= 1;
  }
  return ret;
}

int main(int argc, char *argv[])
{
  int m;
  cin >> m;
  map<long long, int> p;

  for (int i = 0; i < m; i++) {
    long long x;
    cin >> x;
    ++p[x];
  }

  long long ret = 1;
  long long d = 1;
  for (const auto& pr : p) {
    long long x = pr.first;
    long long k = pr.second;

    long long f = mpow(x, k*(k+1)/2, MOD);
    ret = mpow(ret, k+1, MOD) * mpow(f, d, MOD) % MOD;
    d = d * (k+1) % (MOD-1);
  }
  
  cout << ret << endl;
  
  return 0;
}

[Coursera] Game Theory Final

　Courseraのゲーム理論の期末試験を受けた。

授業のクイズ、週ごとの課題がすべて満点だったので、
期末試験も満点取ったるぞーという意気込みで、
土曜日まじめに復習（講義ノートの見直し、課題の見直し）して、日曜テストを受けた。

結果。。。

満点でしたーーヽ( ´￢`)ノ ﾜ～ｲ !!

以下履修した感想を簡単に書いておく。

• 英語は比較的聞きやすかった（字幕なしで理解できた）
• 一部難しい内容があった（Mixed Strategyのナッシュ均衡を求めるアルゴリズムのところ、ベイジアンゲームの定式化）が、概ね易しい内容だった
• 内容を理解するには、大学入試〜大学１・２年くらいの数学力は必要と思われる
• 例題が多く、抽象的な概念を理解するのを助けてくれた
• 課題やクイズは 専門用語の定義を抑えていれば概ね解けると思われる
• 授業の評価とは関係のないゲームがあって、それがなかなか面白かった

モンテカルロ法でシャープレイ値を計算

　昨日動的計画法でシャープレイ値を計算してみたが、一般的な計算方法としてよく知られたものだったらしい。他にもモンテカルロ法を使った計算方法もあるとのことなので、モンテカルロ法で計算するプログラムも書いてみた。

問題は昨日のやつと同じく国連理事国のやつ。random_shuffle便利だよねーっていう。

#include <iostream>
#include <vector>

using namespace std;

const int N = 15;  // # of members
const int M = 9;   // required # of members to approve a decision
const int K = 5;   // # of permanent members
const int T = 1e6; // # of MC iteration

inline int payoff(int s) {
  int vetoMask = (1<<K) - 1;

  if ((s & vetoMask) < vetoMask)
    return 0;

  if (__builtin_popcount(s) < M)
    return 0;

  return 1;
}

double sharpleyValue(int p) {
  vector<int> v;
  for (int i = 0; i < N; i++)
    v.push_back(i);

  double ret = 0;
  for (int _ = 0; _ < T; _++) {
    random_shuffle(v.begin(), v.end());
    int s = 0;
    for (auto & x : v) {
      int t = s | 1 << x;
      if (x == p) {
        ret += payoff(t) - payoff(s);
        break;
      }
      s = t;
    }
  }
  
  return ret / T;
}

int main() {

  for (int i = 0; i < N; i++)
    cout << i << " :" << sharpleyValue(i) << endl;

  return 0;
}

国連理事国のシャープレイ値を計算

　最近シャープレイ値なるものを学んだ。協力ゲームにおいて、ゲームでえられた利得をグループ内のメンバーで公正に分配する方法の一つだ。直感的にはゲームにより貢献したプレイヤーほど多くのシャープレイ値が割り当てられる。

　このシャープレイ値の計算がなかなか面白い。ということでシャープレイ値を計算するプログラムを書いてみた。

[問題設定]

国連の理事国には、常任理事国と非常任理事国が存在する。

常任理事国は５カ国、非常任理事国は１０カ国である。

国連において何かを決める場合、常任理事国すべての賛成に加え、合計９カ国の賛成が必要である。

いま議論している内容が可決となる場合は利得1、否決となる場合は利得0とする。

常任理事国、非常任理事国のシャープレイ値を求めよ。

　以下C++で書いたコード。なんとなくTSPのビットDP解法に似ている気がする。

#include <iostream>

using namespace std;

const int N = 15;  // # of members
const int M = 9;   // required # of members to approve a decision
const int K = 5;   // # of permanent members
long long fac[N+1];

void init() {
  fac[0] = 1;
  for (int i = 1; i <= N; i++)
    fac[i] = i * fac[i-1];
}

inline int payoff(int s) {
  int vetoMask = (1<<K) - 1;

  if ((s & vetoMask) < vetoMask)
    return 0;

  if (__builtin_popcount(s) < M)
    return 0;

  return 1;
}

double sharpleyValue(int p) {

  double ret = 0;
  
  for (int i = 0; i < 1<<N; i++) {
    if (!(i & 1<<p))
      continue;
    
    double tmp = payoff(i) - payoff(i & ~(1<<p));
    int s = __builtin_popcount(i);
    ret += tmp * fac[s-1] * fac[N-s];
  }

  return ret / fac[N];
}

int main() {

  init();

  for (int i = 0; i < N; i++)
    cout << i << " :" << sharpleyValue(i) << endl;

  return 0;
}

結果は、
常任理事国 = 0.1963
非常任理事国　= 0.001865

ということで常任理事国の影響度が非常に高いことが分かる。

[Coursera] Game Theory Week 4: Extensive-Form Games

  Took the 4th week's lecture of "Game Theory." I am half way to the completion and the course is getting difficult now.

Results:

• Video Lectures: All seen and understood
• In-Video Quizz 5/5 (Passed)
• Problem Sets 6/6 (Passed)

What I Learned:

• normal form game vs extensive form game
• perfect information game vs imperfect information game
• pure strategy and Nash equilibrium in extensive form game
• induced normal form 
• sub-game perfect equilibrium
• back induction for finding sub-game perfect equilibrium
• mix strategy vs behavioral strategy

[Coursera] Game Theory Week 3: Alternate Solution Concepts

"N Seoul tower"

 I started taking "Game Theory" lecture of Stanford University on Coursera platform two weeks ago. I am going to keep track of the lectures here on my blog.

Result:

• Video Lectures: all seen and understood
• In-Video Quizzes: 2/2 (passed)
• Problem Sets: 5/5 (passed)

What I learned:

• strictly dominated strategy
• iterated removal of strictly dominated strategies
• relation of Nash Equilibrium and iterated removal method
• maxmin strategy
• minmax strategy
• Minimax theorem (von Neumann, 1928)

整数nが11で割れるかどうか

　整数nが11で割れるかどうかを考える。

10の冪乗の法11での剰余を見てみると

1 = 1 (mod 11)
10 = -1 (mod 11)
100 = 10 * 10 = 1 (mod 11)
1000 = 100 * 10 = -1 (mod 11)
10000 = 1000 * 10 = 1 (mod 11)
.....

のように1,-1,1,-1,...の繰り返しとなる。

よって、nの十進表記を"gfedcba"とするとnの法11での剰余は、

n = g - f + e - d + c - b + a (mod 11)

となる。

つまり、g - f + e - d + c - b + a = 0

であればnは11で割り切れる。

例）
3641 = - 3 + 6 - 4 + 1 = 0 (mod 11)より、3641は11で割り切れる．
13574 = 1 - 3 + 5 - 7 + 4 = 0 (mod 11)より、13574は11で割り切れる．

平方数を4で割った余り

平方数判定をする際にいくらか高速化できるテクニックを知ったのでメモしておく。

nが偶数の平方の場合は、
n = 2k * 2k = 4k^2
と表すことが出来る。よって4で割った余りは0。

nが奇数の平方の場合は、
n = (2k + 1) * (2k + 1) = 4k^2 + 4k + 1
と表すことが出来る。よって4で割った余りは1。

ということで、3とビットのandを取って1より大きくなった場合、その 数は平方数ではないことが分かる。

平方剰余の根

ある整数xに対して

x^2 = a (mod p)

となるようなaを「pの平方剰余」という。

以下ではaが与えられたときに

x^2 = a (mod p)

を満たすxを求めることを考える。ただし今回は簡単のためpが素数の場合のみを考える。

解の存在

ルジャンドル記号(a/p)を以下のように定義する。

a = 0 (mod p)ならば, (a/p) = 0
aが法pで平方剰余ならば、(a/p) = 1
aが法pで平方剰余でなければ、(a/p) = -1

pが奇素数の場合は、オイラーの基準により

(a/p) = a^{(p-1)/2}  (mod p)

が成り立つ。よって繰り返し二乗法を用いることで解の存在を簡単に判定出来る。

long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

解の公式

ルジャンドル記号が1の場合は、

x^2 = a (mod p)

を満たすxが存在する。

p=2の場合はxを全列挙すればよいため、以下ではpは奇素数とする。
pを8で割ったときの余りが3, 5, 7の場合は解の公式が存在する。

i) p = 3 (mod 8)、またはp = 7 (mod 8) の場合

x = ± a^{(p+1)/4} (mod p)

が解となる。

ii) p = 5 (mod 8)の場合

a^{(p-1)/4} = 1なら x = ± a^{(p+3)/8} (mod p)
a^{(p-1)/4} = -1ならx = ± 2^{(p-1)/4} a^{(p+3)/8} (mod p)

が解になる。

// calculate x such that x^2 = a (mod p)
// p is an odd prime that satisfies p = 3, 5, 7 (mod 8)
long long modsqrt(long long a, long long p) {
    
    // case i)
    if (p % 4 == 3)
        return modpow(a, (p+1)/4, p);

    // case ii)
    if (p % 8 == 5) {
        if (modpow(a, (p-1)/4, p) == 1)
            return modpow(a, (p+3)/8, p);
        else
            return modpow(2, (p-1)/4, p) * modpow(a, (p+3)/8, p) % p;
    }
        
    throw;
}

これで75%の場合は公式を使えば解けることが分かる。
しかし残りの25%の場合（p = 1 (mod 8)のとき）は公式が存在しないので、以下で示すアルゴリズムを使って計算するとよい。

非決定性アルゴリズム

Tonelli–Shanks algorithmを使うと、任意の奇素数pについて

x^2 = a (mod p)

を解くことが出来る。

以下に実装例を載せておく。

#include <iostream>
#include <random>

using namespace std;

long long modpow(long long x, long long p, long long mod) {
    long long ret = 1;
    while (p) {
        if (p & 1)
            ret = ret * x % mod;
        x = x * x % mod;
        p >>= 1;
    }
    return ret;
}

// a is a quadratic residue modulo p?
int Legendre(long long a, long long p) {

    long long ret = modpow(a, (p-1)/2, p);

    if (ret == p-1)
        ret = -1;

    return ret;
}

// Tonelli–Shanks algorithm
// calculate R such that R^2 = n (mod p)
// p is an odd prime
// n is a positive integer that satisfies (n/p) = 1
// See http://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm
long long modsqrt(long long n, long long p) {

    // step 1
    long long Q = p-1;
    long long S = 0;
    while (Q % 2 == 0) {
        ++S;
        Q /= 2;
    }

    if (S == 1) {
        return modpow(n, (p+1)/4, p);
    }

    // step 2
    default_random_engine generator;
    uniform_int_distribution<long long> distribution(0, p);
    
    long long z = 1;
    while (Legendre(z, p) != -1) {
        z = distribution(generator);
    }

    long long c = modpow(z, Q, p);
    
    // step 3
    long long R = modpow(n, (Q+1)/2, p);
    long long t = modpow(n, Q, p);
    long long M = S;

    // step 4
    while (t != 1) {
        long long i;
        long long t2 = t;
        for (i = 1; i < M; i++) {
            t2 = t2 * t2 % p;
            if (t2 == 1)
                break;
        }

        long long b = modpow(c, 1LL<<(M-i-1), p);
        R = R * b % p;
        t = (t * b % p) * b % p;
        c = b * b % p;
        M = i;
    }
    
    return R;
}

// test
bool isPrime(long long n) {
    if (n < 2)
        return false;
    
    for (long long i = 2; i * i <= n; i++)
        if (n % i == 0)
            return false;
    
    return true;           
}

#include <cassert>

int main(int argc, char **argv) {
    
    for (int i = 3; i < 10000; i++) {
        if (!isPrime(i))
            continue;

        for (int j = 0; j < i; j++) {
            if (Legendre(j, i) != 1)
                continue;

            long long x = modsqrt(j, i);
            assert (x * x % i == j);
        }

    }

    return 0;
}

参考URL

[1] Euler's criterion - Wikipedia, the free encyclopedia
[2] Computing square roots mod p
[3] Tonelli–Shanks algorithm - Wikipedia, the free encyclopedia

Lucasの定理

　Lucasの定理というものを知った。大きな数の重複なし組み合わせの値を素数を法として求めるときに使える。

非負整数m, n, および素数pについて以下が成り立つ。

 

ただし、m_i, n_iはそれぞれ以下を満たす。

証明はここが分かりやすい。

この定理を見て思ったのは、ほとんどゼロになるのでは？ということだった。
例えば、p=2の場合を考える。m_i < n_i となるようなiが存在すると右辺が0になってしまう。ということは、すべてのiについてm_i ≥ n_iである場合のみ右辺は非零となる。この確率は10桁の二進数だと、pow(3/4, 10) ~ 0.0563程度である。

本当か？ということで試してみた。

#include <iostream>

using namespace std;

int comb[1024][1024];

int main(int argc, char **argv) {
    
    comb[0][0] = 1;
    for (int i = 1; i < 1024; i++) {
        comb[i][0] = comb[i][i] = 1;
        for (int j = 1; j < i; j++) {
            comb[i][j] = (comb[i-1][j-1] + comb[i-1][j]) % 2;
        }
    }

    int total = 0;
    int nonzero = 0;

    for (int i = 0; i < 1024; i++) {
        for (int j = 0; j < 1024; j++) {
            ++total;
            nonzero += comb[i][j];
        }
    }

    cout << 1.*nonzero/total << endl;

    return 0;
}

上のプログラムを実行すると、確かに0.0563と出る。

同様に計算すると、一般にp進数d桁に収まる非負整数m, nでcomb(m, n) % pを計算すると、非零になる割合はpow(1/2 + 1/(2p), d)程度だと考えられる。