n通の便せんとn通の封筒が用意されている。正しい組合せになっていない、つまり、どの便せんも間違った封筒に入れてしまった。そういう可能性は何通りあるか、というのが乱列 (derangement) という。ほかに、例として、預かった荷物がすべて間違ってｎ人に返された場合とか。

n個のものの乱列の数（組合せ数）をd_nと書くと、その計算式は以下の通り与えられている。

n=0からの乱列数は以下の通り。検証に使うといいだろう。

1, 0, 1, 2, 9, 44, 265, 1854, 14833, 133496, 1334961, 14684570, 176214841, 2290792932, 32071101049, 481066515734, 7697064251745, 130850092279664, 2355301661033953

16進数100万桁表示の円周率計算プログラムについて考える。実行速度は2秒以内としたい。効率のよい計算式と、FFTによる乗算で高速化にチャンレンジしたい。

ここのブログでも紹介している多倍長整数を利用するが、16進数に合わせて、1語の大きさを2進32ビット、つまり、16進8桁に直す。10進数と違って、除算ではビットシフト、剰余計算ではビットマスクだけでできるので、さらなる高速化が期待できそう。

整数一つの配列はこんなものになるだろう。オーバーフローが起きて、64ビット整数にしないといけない可能性もあるが。

　　unsigned num[125000];

Chudnovsky)の公式は計算が速いようだが、確かめないとなんともいえない。

ネット上を調べたら、全ての桁をいっぺんに計算するのではなく、特定の桁だけを算出する方法もあるらしい。今回の用途に合致するかもしれない。

英文 Algorithm for Calculating Individual Hexadecimal Digits of Pi

なお、16進円周率の値一部はここにて確認できる。

オイラーの関数φ(n)の値が与えられた時の、nを求める問題。つまり、

　 φ(n)=x のxに対し、n=φ^-1(x)

を考える。

整数の素因数分解と同様、逆関数の計算はとても難しい。だから問題として出されているかもしれない。

1. 元問題 — UVa #11073 – Euler’s Totient Function
　　xの値の範囲は 1<=x<=1000000000。
nの値が複数ある場合、昇順で表示する必要もある。

2. 実例
　　x=2に対して、n=3,4,6
　　x=4に対して、n=5,8,10,12
　　x=583925760に対して、書ききれないが、異なるnは19639個存在する。

3. 考え方
　xの値が小さければ、例えば、100万までなら、φ(n)の値をすべてテーブルに計算登録し、xでlookupする方法も考えられるが、10億もの巨大xではこの方法は無理だろう。

幸い、元問題の入力値xの個数はそれほど多くないらしく、素直に個々にxの逆関数を計算していくことにした。

自分では以下の漸近式を使った。

つまり、2から約32000（10億の平方根）までのの素数をひとつひとつ持ってきて、xを(p-1)*p^kで割り切れるかどうかを確かめる。割り切れるのであれば、x/((p-1)*p^k)を新たなxとして、その逆関数を再呼び出して分解していく。計算できたnをひとつひとつテーブルに登録し、最後にソートをかけ画面出力する。

ちなみに、本方法では、途中で打ち切ることがあっても、nの値がダブって算出されることはない。

なお、x=1に対しては、n=1,2と結果表示し、他の奇数のxに対しては、解がなしと表示する。

4. 実例で説明
　基本的な考え方は上記の通りだが、細かいところについては以下の実例をよく見て、理解すること。

x=4について考える。

(1-1) まず最初の素数p=2を持ってくる。k=0とし、(p-1)*p^kを見ていく。なお、k=1については(1-2)の所、k=2については(1-3)の所で調べる。
　　p=2, k=0では、(p-1)*p^k=1, p^k+1=2, x/((p-1)*p^k))=4になり、φ^-1(4)→2*φ^-1(4)になる。
　(1-1-1) 素数2は上で既に使ったので、ここでは次の素数3で調べる。p=3,k=0のもとでは、(p-1)*p^k=2, p^k+1=3, 4/((p-1)*(p^k)=2になり、2*φ^-1(4)→6*φ^-1(2)になる。
　　(1-1-1-1) つぎの素数5で調べる番になっているが、2を(p-1)*p^k=4*5^kで割り切ることはありえないので、一旦ここで打ち切り。
　(1-1-2) つぎの素数5で調べる。p=5,k=0のもとでは、(p-1)*p^k=4, p^k+1=5, 4/((p-1)*(p^k)=1になり、2*φ^-1(4)→10*φ^-1(1)になる。φ^-1(1)は1であるので、最初のnが見つかり、n=10。
(1-2) (1-1)の続きをやる。つまり、p=2, k=1。
　　p=2, k=1では、(p-1)*p^k=2, p^k+1=4, x/((p-1)*p^k))=2になり、
φ^-1(4)→4*φ^-1(3)になる。
　(1-2-1) 次の素数3で調べる。p=3,k=0のもとでは、(p-1)*p^k=2, p^k+1=3, 4/((p-1)*(p^k)=2になり、
4*φ^-1(2)→12*φ^-1(1)になり、2番目のn=12がみつかる。
(1-3) (1-2)の続きをやる。つまり、p=2, k=2。
　　p=2, k=2では、(p-1)*p^k=4, p^k+1=8, x/((p-1)*p^k))=1になり、
φ^-1(4)→8*φ^-1(1)になり、3番目のn=8がみつかる。

(2-1) 2番目の素数p=3を持って来て調べる。
　　p=3, k=0では、(p-1)*p^k=2, p^k+1=3, x/((p-1)*p^k))=2になり、
φ^-1(4)→3*φ^-1(2)になるが、つぎの素数5ではそれ以上nを見つけることはないので、打ち切る。

(3-1) 3番目の素数p=5を持って来て調べる。
　　p=5, k=0では、(p-1)*p^k=4, p^k+1=5, x/((p-1)*p^k))=1になり、
φ^-1(4)→5*φ^-1(1)になり、4番目のn=5がみつかる。

ということで、なかなか理解するのは大変かもしれないが、個々の素数をつぎつぎと持ってきてテストし、つぎの深いステップでは、同じ素数は使ってはいけないところが味噌。また、同じステップのところでは、kの値も0からインクリメントして調べる。いわば、3次元的に調べるイメージになる。

5. 結果
あまりよい計算方法がないらしく、素数をひとつひとつで確かめていき、効率の面を心配していた。結果的にランク1位の実行速度になり、皆も相当苦労していたのだなと実感した。

まだまだ研究余地の残っている問題だと思う。

オイラーの関数φ(n)とは、1からnまでの整数のうちnと互いに素なものの個数のこと。

たとえば、1,2,3,4,5,6のうち6と互いに素なのは、1,5の二つなので、φ(6)=2。素数pについては、1からpまでの整数のうち互いに素でないのはp自身だけなので、φ(p)=p-1となることは直ちに理解できよう。なお、φ関数は英語では、totient functionと呼ばれている。

1からのφ関数の値は以下の通り。

　1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12

φ(n)の計算に重要なもの（性質）をいくつかピックアップしておく。

整数の素因数分解と同様、個々のφ(n)の計算と、多数のnのφ(n)の計算とを分けて、その計算法を考えよう（数学ではそんなことを考える必要性は全くないが）。

個々のφ(n)の計算では、nの素因数分解を利用する方法が最も効率的だろう。nに含まれている素因数p_kが分かれば、(n/P_k)*(p_k-1)で計算していけばよし。

多数のnをまとめてそのφ(n)を計算するなら、以下のプログラムを利用しよう。φ(n)の性質（互いに素）を利用した方法。

#define MAX 1000000
int phi[MAX+10];
char prime[MAX+10];
void euler_phi(void)
{
    register int i, j;
    for (i = 0; i < MAX; i++) phi[i] = i;
    phi[2] = 1;
#if 0
    memset(prime, 0, sizeof(prime));
#endif
    for (j = 4; j < MAX; j += 2) prime[j] = 1, phi[j] -= phi[j] >> 1;
    for (i = 3; i < MAX; i++) {
        if (prime[i]) continue;
        phi[i] -= phi[i] / i;
        for (j = i+i; j < MAX; j += i) prime[j] = 1, phi[j] -= phi[j] / i;
    }
#if 0
    for (i = 1; i <= 100; i++) printf("(-452953008,-2082935264) ", i, phi[i]);
    printf("\n");
#endif
}

これでいっぺんに、MAX（100万）までの値が計算される。

LIS (Longest Increasing Subsequence) とは最長単調増加シーケンスのこと。数字列（あるいは文字列等、互いに大きさの比較ができるものなら何でもOK）の中から、値が大きくなっていき、もっとも長いシーケンス（数字列）を見つけ出すアルゴリズムのことです．

　-7, 10, 9, 2, 3, 8, 1

例えば，上の数字列では，LIS は -7 2 3 8 になり、その長さは4だ。

また、概念は同じだが、LDS (Longest Decreasing Subsequence、最長単調減少シーケンス) というのもある。

ここでは、i (1<=i<=N、Nはシーケンスの長さ）番目から始まるLIS（またはLDS）の長さについて考える。

上の例では、LIS[i] = { 4, 1, 1, 3, 2, 1, 1 }
　　　　　　LDS[i] = { 1, 3, 2, 2, 2, 1, 1 } となる。

各データの位置から、LIS(LDS)の長さが分かれば、LISシーケンスを作り出すことができる。

#define SIZE 10000 /* 必要に応じて変更 */
int N; /* データの個数 */
int num[SIZE+3]; /* データ */
int LIS[SIZE+3]; /* LIS長さ */
int LDS[SIZE+3]; /* LDS長さ */
void calc(void)
{
    int i, j;
    int mi, md;
    LIS[N-1] = LDS[N-1] = 1;
    for (i = N-2; i >= 0; i--) {
        mi = md = 0;
        for (j = i+1; j < N; j++) {
            if (num[i] < num[j]) {
                if (LIS[j] > mi) mi = LIS[j];
            } else {
                if (LDS[j] > md) md = LDS[j];
            }
        }
        LIS[i] = mi+1;
        LDS[i] = md+1;
    }
}

計算時間は O(N^2) 。

下のコードはO(nlogn)の高速バージョン。LISのみになっているが、LDSへの対応もそれほど困難ではないだろう。

#define SIZE 10005  /* 配列のサイズ。適宜変更 */
int N;  /* データの個数 */
int max;  /* 計算終了の時点では最長LISの値が入る */
int num[SIZE+3];  /* データ */
int LIS[SIZE+3];  /* それぞれのデータ置からのLIS長さ */
int tmp[SIZE+3];  /* 作業用 */
void calc(void)
{
    int k;
    int lo, hi, mid;
    LIS[N-1] = 1;
    tmp[0] = num[N-1];
    max = 1;
    for (k = N-2; k >= 0; k--) {
        if (tmp[max-1] > num[k]) { tmp[max++] = num[k]; LIS[k] = max; continue; }
        if (tmp[max-1] == num[k]) { LIS[k] = max; continue; }
        lo = 0;
        hi = max;
        while (lo < hi) {
            mid = (lo + hi) >> 1;
            if (tmp[mid] > num[k]) lo = mid+1;
            else hi = mid;
        }
#if 0
        if (lo < max) tmp[lo] = num[k], LIS[k] = lo+1;
        else tmp[max++] = num[k], LIS[k] = max;
#else
        tmp[lo] = num[k], LIS[k] = lo+1;
#endif
    }
}

問題意識は UVa #11480 Jimmy’s Balls による。

つまり、自然数Nを３つの自然数 x, y, z に分割する時の分割数（組合せ数）。ただし、x < y < z。自然数とは1以上の整数のこと。

たとえば、N=10に対して、(x, y, z)は、(1,2,7), (1,3,6), (1,4,5), (2,3,5) の４通りしか選べないので、分割数は4になる。

割とよく遭遇しそうな問題だが、一般解は以下の数式で表される。

　　分割数 f = 3p(p-1)+1　　if m = 0

　　　　　 = 3p(p-1)+mp　　if m > 0

ただし、p=N/6 （6で割った時の商）、m=N%6（6で割った時の余り）。

いくつかの検証データ(N,f)を書いておく。

　(5,0) (6,1) (7,1) (8,2) (9,3) (10,4) (11,5) (12,7) (13,8) (14,10) (15,12) (100,784) (500,20584)

一般化して、Nをn個に分割する時の数式ってあるのだろうか。