平方数について

平方数について簡単なまとめ。

1.平方数とは？

整数の二乗で表すことのできる数のこと。0, 1, 4, 9, 16, ...と続く。

2. ２つの平方数の和

２つの平方数の和に関しては、フェルマーのクリスマス定理が有名。この定理では奇素数が２つの平方数の和で表すことができるか否かについて述べられているが、

(x₁² + y₁²)(x₂² + y₂²) = (x₁x₂ - y₁y₂)² + (x₁y₂ + y₁x₂)²

の公式より、一般の自然数nについて以下のことが言える。

n = n₁² * n_2, ただし、n₂は法4の下で3と合同な奇素数を因数にもたない。

と分解できる場合、nは２つの平方数の和で表すことができる。上記のように分解できない場合は２つの平方和で表すことはできない。

3. ３つの平方数の和

３つの平方和ですべての自然数を表すことができるか？答えはできない。法8の下で7と合同な自然数は３つの平方数の和で表すことはできない。なぜなら、平方数は8を法として{0, 1, 4}のいずれかと合同になるが、これらの任意の３つの組み合わせの和は7と合同になりえないからである。

4. ４つの平方数の和

4つの平方数の和を用いれば、2. および3.で表すことができなかった自然数を表すことができる。

7 ≡ 3 (mod 4)および7 ≡ 7 (mod 8)より、7は２つの平方数の和でも３つの平方数の和でも表すことができない。しかし、

7 = 1+1+1+4

と4つの平方数の和で表すことができる。

すべての自然数は、４つ以下の平方数の和で表すことができる。定理の名前は知らなかったが、「ラグランジュの四平方定理」というらしい。

SRM 527 DIV1 450 P8XMatrixRecovery

2部グラフのマッチングまでは、気付いてた。あと一歩だった。
"ある条件を満たす辞書順最小のものを求めよ。"という問題はgreedyで解ける場合が多いような気がする。

view plainprint?

1. int r;  
2. int c;  
3. bool edge[60][60];  
4. bool used[60];  
5. int pr[60];  
6.   
7. class P8XMatrixRecovery {  
8.     bool dfs(int s) {  
9.         used[s] = 1;  
10.   
11.         REP(v, 2*c) {  
12.             if (!edge[s][v])  
13.                 continue;  
14.             int w = pr[v];  
15.             if (w == -1 || (!used[w] && dfs(w))) {  
16.                 pr[s] = v;  
17.                 pr[v] = s;  
18.                 return true;  
19.             }  
20.         }  
21.         return false;  
22.     }  
23.   
24.     int match(vector<string> &rows, vector<string> &cols) {  
25.         memset(edge, 0, sizeof(edge));  
26.         REP(i, c) REP(j, c) {  
27.             edge[i][j+c] = edge[j+c][i] = 1;  
28.             REP(k, r) {  
29.                 if (cols[j][k] != '?' && rows[k][i] != '?' && rows[k][i] != cols[j][k]) {  
30.                     edge[i][j+c] = edge[j+c][i] = 0;  
31.                     break;  
32.                 }  
33.             }  
34.         }  
35.   
36.         int ret = 0;  
37.         memset(pr, -1, sizeof(pr));  
38.         REP(i, 2*c) {  
39.             if (pr[i] >= 0)  
40.                 continue;  
41.             memset(used, 0, sizeof(used));  
42.             if (dfs(i))  
43.                 ++ret;  
44.         }  
45.   
46.         return ret;  
47.     }  
48.   
49. public:  
50.     vector<string> solve(vector<string> rows, vector<string> columns) {  
51.         r = rows.size();  
52.         c = rows[0].length();  
53.   
54.         REP(i, r) REP(j, c) {  
55.             if (rows[i][j] == '?') {  
56.                 rows[i][j] = '0';  
57.   
58.                 if (match(rows, columns) < c)  
59.                     rows[i][j] = '1';  
60.             }  
61.         }  
62.   
63.         return rows;  
64.     }  
65. };  

int r;
int c;
bool edge[60][60];
bool used[60];
int pr[60];

class P8XMatrixRecovery {
    bool dfs(int s) {
        used[s] = 1;

        REP(v, 2*c) {
            if (!edge[s][v])
                continue;
            int w = pr[v];
            if (w == -1 || (!used[w] && dfs(w))) {
                pr[s] = v;
                pr[v] = s;
                return true;
            }
        }
        return false;
    }

    int match(vector<string> &rows, vector<string> &cols) {
        memset(edge, 0, sizeof(edge));
        REP(i, c) REP(j, c) {
            edge[i][j+c] = edge[j+c][i] = 1;
            REP(k, r) {
                if (cols[j][k] != '?' && rows[k][i] != '?' && rows[k][i] != cols[j][k]) {
                    edge[i][j+c] = edge[j+c][i] = 0;
                    break;
                }
            }
        }

        int ret = 0;
        memset(pr, -1, sizeof(pr));
        REP(i, 2*c) {
            if (pr[i] >= 0)
                continue;
            memset(used, 0, sizeof(used));
            if (dfs(i))
                ++ret;
        }

        return ret;
    }

public:
    vector<string> solve(vector<string> rows, vector<string> columns) {
        r = rows.size();
        c = rows[0].length();

        REP(i, r) REP(j, c) {
            if (rows[i][j] == '?') {
                rows[i][j] = '0';

                if (match(rows, columns) < c)
                    rows[i][j] = '1';
            }
        }

        return rows;
    }
};

メルセンヌ素数と完全数の関係

メルセンヌ素数と完全数には面白い関係がある。

1. メルセンヌ素数とは？

　メルセンヌ素数とは、メルセンヌ数かつ素数であるような数のこと。メルセンヌ数とは、2^k-1という形で表される数のことで、2進数で表したときに、1のみで構成される数のことである。

2. 完全数とは？

完全数とは、その数のproper divisor（自分以外の約数）の和が自分自身と等しいような数のこと。

6 = 1 + 2+ 3

28 = 1 + 2 + 4 + 7 + 14

が最初の2つの完全数である。

3. 面白い関係って？

　メルセンヌ素数と完全数の間には以下の定理が成り立つ。

メルセンヌ数2^k-1が素数であるとき、2^k-1(2^k-1)は完全数である。逆に、すべての偶数の完全数は、2^k-1(2^k-1)という形（2の冪乗とメルセンヌ素数の積）で表現できる。（※注意）

　証明のために関数σを導入します。σ(n)は、nの約数（properな約数ではなく、自分自身も含むことに注意）の和を表します。例えば、

σ(6) = 1 + 2 + 3 + 6 = 12

です。nが完全数の場合は、σ(n) = 2nです。

　

　この関数σについて以下の補題を考えます。

N = ab（a, bは互いに素）と表せるとき、σ(N) = σ(a) * σ(b)である。

　証明)　σ(a) * σ(b) = (a1 + a2 + a3 + ... ) * (b1 + b2 + b3 + ... ) = Σ(aの約数 * bの約数)である。ここでa1,a2,..はaの約数、b1,b2,..はbの約数である。よって、（aの約数 * bの約数）がNの約数となることと、Nの任意の約数dがある一組の(i, j)についてd = ai * bjとなることを示せばよい。

　まず前半だが、a, bが互いに素であるため、a x + b y = 1となるような整数(x, y)が存在する。両辺にNをかけると、a x N + b y N = Nとなる。aの約数 | a かつ bの約数 | Nなので、（aの約数 * bの約数）はa x Nを割りきる。同様に（aの約数 * bの約数）はb y Nを割りきる。よって（aの約数 * bの約数）はNを割りきる。

　次に後半を示す。a, bの素因数分解を

a = p₁^k1p₂^k2p₃^k3...p_n^kn

b = q₁^t1q₂^t2q₃^t3...q_m^tm

と表記する。ki > 0, i = 1,2,...n、およびtj > 0, j = 1,2,..mである。また、aとbが互いに素なので、p_i≒q_j, i=1,2,..n, j =1,2,...mである。

N= a * bなので、Nの素因数分解は

N = p₁^k1p₂^k2p₃^k3...p_n^kn q₁^t1q₂^t2q₃^t3...q_m^tm

となる。ここで、Nの任意の約数は、

x = p₁^k'1p₂^k'2p₃^k'3...p_n^k'n q₁^t'1q₂^t'2q₃^t'3...q_m^t'm

と表すことができる。ただし、0 ≦k'i ≦ki, i = 1,2,...n、0 ≦ t'j ≦tj, j=1,2,...,m.である。

Nの任意の約数を選ぶことと、(k'1,k'2, ..., k'n, t'1,t'2,..., t'm)の値を決めることは一対一に対応している。(k'1,k'2, ..., k'n, t'1,t'2,..., t'm)が決まれば、(k'1,k'2, ..., k'n)、(t'1,t'2,..., t'm)が決まり、これは、対応するaの約数、bの約数が一意に決まることを意味する。□

　それでは、メルセンヌ素数と完全数に関する定理の証明に入ります。まず、最初の部分を証明します。

2^k-1は素数なので、σ(2^k-1) = 1 + 2^k-1 = 2^kです。2^k-1は、2のべき乗なので約数は1, 2, 2², .... 2^k-2, 2^k-1である。よって、σ(2^k-1) = 1 + 2 + 2² + ... + 2^k-1 = 2^k-1。2^k-1（2の冪乗）と2^k-1（奇数）は互いに素なので補題より、σ(2^k-1(2^k-1)) = σ(2^k-1) * σ(2^k-1) = (2^k-1) * 2^k = 2 * 2^k-1(2^k-1)となる。

σ(n) = 2nなので、2^k-1(2^k-1)は完全数である。

　次に後半を示す。偶数の完全数を2^k-1mと表記する。ただし、k≧2、mは奇数。

補題よりσ(2^k-1m) = σ(2^k-1) * σ(m) = (2^k-1) σ(m)。

また、2^k-1mは完全数であることから、σ(2^k-1m) = 2 * 2^k-1m = 2^km。

よって、

2^km = (2^k-1) σ(m)。

ここで、(2^k-1) | 2^km、かつ2^kと(2^k-1)は互いに素であることから、(2^k-1) | m。

よって、m = (2^k-1) Mとして先ほどの式に代入すると

2^k (2^k-1) M = (2^k-1) σ(m)となり、両辺を (2^k-1) で割ると、

2^kM = σ(m)。

m + M = (2^k-1) M + M = 2^kMより、σ(m) = m + Mとなるが、これはmが素数で、M=1であることを意味する。よって、m = 2^k-1（素数）となり、もともとの完全数は、2^k-1(2^k-1)と表される。□

（※注意）奇数の完全数が存在するか否かは不明です。現時点で発見されている完全数はすべて偶数です。

参考：

[1] http://primes.utm.edu/notes/proofs/EvenPerfect.html
[2] http://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%AB%E3%82%BB%E3%83%B3%E3%83%8C%E6%95%B0

n進数表記の数を(n-1)で割れるかどうか

面白い問題があったので、紹介。

http://poj.org/problem?id=1152

　n進数表記の数値を(n-1)で割り切れるかどうかは、すべての桁の数値を足した数が(n-1)で割れるかどうかを調べればいいです。例えば、10進数表記の111111111、45、1233、123453は桁の合計値が9で割れるので9の倍数です。

　n進数表記のm桁の数 A = A_(m-1).. A_2 A_1 A_0について考えます。つまり、

A = A_(m-1) * n^(m-1) + ... A_2 * n^(2) + A_1 * n^(1) + A_0

です。 n = 1 (mod (n-1))なので、

A = A_(m-1) + ... + A_2 + A_1 + A_0 (mod (n-1))

です。

ということで、Aの各桁の合計値を出してそれが(n-1)で割れるかどうかをみればOKです。

　C++のスパゲッティなソースを張っておきます。（dividable()の中で毎回合計値計算してるのは無駄ですね。。最初まじめにHormerの公式を使って余りを計算して解こうとしていたのの名残です。）

view plainprint?

1. int num[128];  
2.   
3. bool dividable(int k, const string &str) {  
4.     int len = str.length();  
5.     int digsum = 0;  
6.     REP(i, len) digsum += num[(int)str[i]];  
7.   
8.     return digsum % (k-1) == 0;  
9. }  
10.   
11. int main() {  
12.     string str;  
13.   
14.     for (int i = '0'; i <= '9'; i++) num[i] = i - '0';  
15.     for (int i = 'A'; i <= 'Z'; i++) num[i] = i - 'A' + 10;  
16.     for (int i = 'a'; i <= 'z'; i++) num[i] = i - 'a' + 36;  
17.   
18.     while (cin >> str) {  
19.         int len = str.length();  
20.         int base = 0;  
21.         REP(i, len) base = max(base, num[(int)str[i]]);  
22.         base = max(2, base+1);  
23.   
24.         int k;  
25.         for (k = base; k <= 62; k++)  
26.             if (dividable(k, str))  
27.                 break;  
28.   
29.         if (k <= 62) cout << k << endl;  
30.         else cout << "such number is impossible!" << endl;  
31.     }  
32.   
33.     return 0;  
34. }  

int num[128];

bool dividable(int k, const string &str) {
    int len = str.length();
    int digsum = 0;
    REP(i, len) digsum += num[(int)str[i]];

    return digsum % (k-1) == 0;
}

int main() {
    string str;

    for (int i = '0'; i <= '9'; i++) num[i] = i - '0';
    for (int i = 'A'; i <= 'Z'; i++) num[i] = i - 'A' + 10;
    for (int i = 'a'; i <= 'z'; i++) num[i] = i - 'a' + 36;

    while (cin >> str) {
        int len = str.length();
        int base = 0;
        REP(i, len) base = max(base, num[(int)str[i]]);
        base = max(2, base+1);

        int k;
        for (k = base; k <= 62; k++)
            if (dividable(k, str))
                break;

        if (k <= 62) cout << k << endl;
        else cout << "such number is impossible!" << endl;
    }

    return 0;
}

最大クリーク問題

http://poj.org/problem?id=3692

を解いていて、「最大クリーク問題」という問題について知った。

　与えられたグラフの部分グラフのうち完全グラフとなるものをクリーク（clique）と呼ぶ。節点数が最大となるクリークを求める問題を「最大クリーク問題」という。この問題はNP完全らしい。

　完全グラフの補グラフは、null graph（枝集合が空集合であるグラフ）であるので、あるグラフの最大クリークを求めることと、その補グラフの最大独立集合を求めることは同値である。

　二部グラフにおいては最小点カバーと最大マッチングが一致すること、および、一般のグラフにおいて最小点カバーと最大独立集合の和が節点数と一致することから、対象の補グラフが二部グラフである場合には最大クリークを簡単に求めることができる。

　”クリーク”という言葉は、社会学から生まれたらしい。ソーシャルグラフの中からお互いが全員知り合いであるような最大集合の大きさを求めたいというのが研究のはじまり。まさに、PKUの問題もこれと同様の問題。実際のソーシャルグラフの補グラフは必ずしも二部グラフになるとは限らないため、さまざまな近似アルゴリズムが研究されているらしい。

　PKUの問題の解答例を載せておきます。

view plainprint?

1. int n;  
2. bool edge[400][400];  
3. bool used[400];  
4. int pr[400];  
5.   
6. bool match(int s) {  
7.     used[s] = true;  
8.     REP(i, n) {  
9.         if (!edge[s][i]) continue;  
10.         if (pr[i] == -1 || (!used[pr[i]] && match(pr[i]))) {  
11.             pr[s] = i;  
12.             pr[i] = s;  
13.             return true;  
14.         }  
15.     }  
16.   
17.     return false;  
18. }  
19.   
20. int main() {  
21.     int b, g, m;  
22.   
23.     int t = 0;  
24.     while (cin >> b >> g >> m) {  
25.         if (!b && !g && !m) break;  
26.   
27.         int x, y;  
28.         n = b + g;  
29.         REP(i, n) REP(j, n) edge[i][j] = 1;  
30.         REP(i, b) REP(j, b) edge[i][j] = 0;  
31.         REP(i, g) REP(j, g) edge[i+b][j+b] = 0;  
32.   
33.         REP(i, m) {  
34.             cin >> x >> y;  
35.             --x, --y;  
36.             edge[x][y+b] = edge[y+b][x] = 0;  
37.         }  
38.   
39.         int ret = n;  
40.         memset(pr, -1, sizeof(pr));  
41.         REP(i, n) {  
42.             if (pr[i] == -1) {  
43.                 memset(used, 0, sizeof(used));  
44.                 if (match(i)) --ret;  
45.             }  
46.         }  
47.   
48.         printf("Case %d: %d\n", ++t, ret);  
49.     }  
50.   
51.     return 0;  
52. }  

int n;
bool edge[400][400];
bool used[400];
int pr[400];

bool match(int s) {
    used[s] = true;
    REP(i, n) {
        if (!edge[s][i]) continue;
        if (pr[i] == -1 || (!used[pr[i]] && match(pr[i]))) {
            pr[s] = i;
            pr[i] = s;
            return true;
        }
    }

    return false;
}

int main() {
    int b, g, m;

    int t = 0;
    while (cin >> b >> g >> m) {
        if (!b && !g && !m) break;

        int x, y;
        n = b + g;
        REP(i, n) REP(j, n) edge[i][j] = 1;
        REP(i, b) REP(j, b) edge[i][j] = 0;
        REP(i, g) REP(j, g) edge[i+b][j+b] = 0;

        REP(i, m) {
            cin >> x >> y;
            --x, --y;
            edge[x][y+b] = edge[y+b][x] = 0;
        }

        int ret = n;
        memset(pr, -1, sizeof(pr));
        REP(i, n) {
            if (pr[i] == -1) {
                memset(used, 0, sizeof(used));
                if (match(i)) --ret;
            }
        }

        printf("Case %d: %d\n", ++t, ret);
    }

    return 0;
}

参考：
1. simezi_tanの日記
2. Clique problem - wikipedia

サーバー奮闘記（18） gitを使う

"Cowboy Cafe - Dubois - Wyoming - July 2005" By GeordieMac Pics

　結構前にgitはローカルとvpsともに入れていたのですが、GitHubにソースを上げるときくらいしか使ってませんでした。

　vpsにマスターのレポジトリを作って、ローカルにいくつかレポジトリを作ってpushしたりpullしたりして遊んでみました。簡単なメモを残します。

1. sshのconfigファイルの設定

　sshを使ってリモートログインしたいサーバーが複数台ある場合は、configファイルを書いておくと便利です。サーバーのアドレスや、sshのポート番号、どのrsa秘密鍵を使うのか、などを記述しておくと、sshコマンドを使う際にいろいろとオプションを書かなくてよくて楽です。

~/.ssh/config

HOST vps

　HostName 123.12.1.123

　Port 4321

　IdentityFile ~/.ssh/id_rsa_vps

　User taro

　それぞれのサーバーに付き、上のような設定を書いておくと、

ssh vps

のようにHOST名を書くだけでログインできるようになります。

2. サーバー側にbareリポジトリを作る

　調べてみた感じだと、2つやり方があるようです。1つ目は、サーバー上で普通のレポジトリを作成して、それをbareリポジトリにcloneして、さらにそれをローカルにcloneするという方法。2つ目は、サーバー側にbareレポジトリを作成して、クライアント側でリモートサーバーを登録する方法です。

1つ目の方法だと、cloneした時点で、push/pull先が自動でremoteに登録されるので明示的にリモートサーバー名を追加する必要はありません。ただし、サーバー側でcloneする前に、適当なファイルを作ってcommitしておかないといけません。

参考: http://d.hatena.ne.jp/Naruhodius/20110415/1302832214

2つ目の方法は、サーバー側では

$ git --bare init

とするだけです。ただし、1つ目の方法と違って、remoteサーバー名が自動で登録されないので、サーバーに名前を付けたい場合（※）は、

$ git remote add origin ssh://vps/~taro/git/repos/hello

のようにremote サーバーを明示的に追加する必要があります。

参考: http://blog.champierre.com/670

※remoteサーバー名は必ず登録する必要はないですが、毎回push/pullする際に

$ git push ssh://vps/~taro/git/repos/hello

とするのは面倒なので、登録しておいた方が便利。

おそらく2つ目のやり方の方が一般的かと思います。