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

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

derange.jpg

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

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

080907.jpg

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

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

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

  unsigned num[125000];

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

p11302.jpg

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

英文 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の逆関数を計算していくことにした。

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

p11073.jpg

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

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

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

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

x=4について考える。

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

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

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

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

p11073-1.jpg

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)の計算に重要なもの(性質)をいくつかピックアップしておく。

totient.jpg

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

個々のφ(n)の計算では、nの素因数分解を利用する方法が最も効率的だろう。nに含まれている素因数pkが分かれば、(n/Pk)*(pk-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万)までの値が計算される。

等差数列に関わる問題。

1. 問題文
長さp (3<=p<=30内の素数)の等差数列にならない最小数列(細かい定義は元問題文を参照)に対し、そのn (1<=n<=20000000000) 番目の項を示せ。

2. 実例
p=3の時の最小数列は 1,2,4,5,10,11,13,14,28,29,… になっている。数列のどの3項を取ってきても等差数列にならない。

しかし、たとえば、1,2,4,5,9,10,… という数列ではダメだ。1,5,9が等差数列になってしまうから。

3. 解き方
ブログ内の昔の記事にも書いたが、数列を扱う問題では、オンライン整数列大辞典が大いに参考になる。

本問題に関しても、1,2,4,5,10,11と数列を打ち込むと直ちに ID: A003278が表示される。COMMENTを読み、とくにFORMULAを注意深く観察すると、それらしき計算式が表示されていることが多い。今回も

a(n) = b(n+1) with b(0)=1, b(2n)=3b(n)-2, b(2n+1)=3b(n)-1

という一行が役に立ちそう。

数式に沿ってプログラムを組み、最初に数項が正しく計算されることを確認し、3以外の素数への拡張を推測していく。省略するが、最終的に以下の通り計算を行う。

int p;  /* 問題文の素数 */
long long calc(long long n)
{
    if (n == 0) return 1;
    return p*calc(n/(p-1)) + n%(p-1) - (p-1);
}

上の関数を calc(n-1) で呼び出す。nの値の対数回計算で結果が出る。

問題文を楽しむなら、ここで紹介したやり方ではいけないが、自分で同じ数式を探し当てることができるなら素晴らしいけど、無駄なことをやらなくてもいいじゃないかな。

計算速度が至上命令か 。本問題のプログラム実行制限時間はついに0.666秒に設定されてしまった。自分としては嫌なんだね、はっきりと。

アルゴリズムの優劣性については、実行速度がとても大事だが、可読性、そこからの発展性、理解しやすさ、など場面や用途によって様々なものがあっていいはず。ソートアルゴリズムについても、O(N)のものがあっても、ソートの本質を理解する上では、様々なものが用意されている。実行制限時間をきつく設定したら、アルゴリズムの強制採用にほかならない。大反対だね。

1. 問題文
n (1<=n<=10000) 個のデータ (整数、10進8桁以内、負数もありえる)が与えられている。さらに、様々な値m (1<=m<=n) が与えられ、m個のLISシーケンス(Longest Increasing Subsequence、最長単調増加シーケンス)をデータからピックアップして表示せよ、というのが問題の意味。なお、LISシーケンスはなるべくデータの先頭にあるものを選択しないといけない。

2. 実例
例えば、以下のn=10個のデータがあったとする。m=1,2,3,…10のLISシーケンスはそれぞれ以下の通り。

-3 19 12 -7 18 23 -5 0 11 22

m=1 LIS=-3

m=2 LIS=-3 19

m=3 LIS=-3 19 23

m=4 LIS=-3 12 18 23

m=5 LIS=-7 -5 0 11 22

m>5 LISは存在しない

3. 計算法
O(nlogn)バージョンのLISアルゴリズムを使う。ブログに紹介があるので、参照されたい。

つまり、n個のデータに対するそれぞれのLIS長をまず算出する。上の実例では、次になる。

data (-3 19 12 -7 18 23 -5 0 11 22)

LIS (4 2 3 5 2 1 4 3 2 1)

最長LISは5という結果も上のLISテーブルから分かる。

つぎに、それぞれのmに対して、そのLISシーケンスをピックアップする。

m=1 LIS1=4 >= 1 なので、data1=-3がその結果。
m=2 LIS1=4 >= 2、LIS2=2 >= 1 しかも data2=19 > data1なので、-3 19が結果になる。
m=3 LIS1=4 >=3、LIS2=2 >=2 しかも data2 > data1、data2=19がLISシーケンス2番目のものになる。さらに、LIS3=3 は 1以上だが、値はdata2以下なので、却下。LIS4, 5も同様の理由で却下。LIS6=1が1以上だし、値data6=23は19よりも大きいので、3番目はdata6になる。
m=4については、data2=19のLIS2は2なので、3よりも小さいことから、2番目にはならない。隣のdata3=12が2番目になる。つぎの12よりも大きく、LISの値が2以上のものは18なので、3番目になる。4番目は18よりも大きく、LISの値が1以上のもの、つまり、23となる。
m=5については省略するが、先頭からまずLISの値が5よりも大きいものを探し、-7を取る。つぎに-7よりも値が大きく、LISの値が4よりも大きいものを取る。

全体の実行時間は、LISテーブルの算出にO(nlogn)、LISシーケンスの抽出にO(n)でできるので、全体でもO(nlogn)になる。

4. 結果
入力部分を工夫して、最終的にランク1位の実行速度と同じにした。出力部分はprintf()をつかったままなので、まだ余力は残っているが、誤差範囲内の時間(0.01秒)を縮めても意味ないので、そのままにしておいた。

ともかく、スピードが命という風潮に強く警告したい気持ちでいっぱい。いくらなんでもやり過ぎだな。

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

p11271.jpg

XY平面に無限に広がる抵抗ネットワークに対し、任意の両点(格子点)間の抵抗値を求める問題。自分は大学では電子工学を勉強したのに解けない。情けないことかな。

1. 問題
任意の両点だが、無限に広がっているので、原点をどこに取っても同じはず。そこで、両点の片方を原点にし、原点からの抵抗値を求めるというのが本問題。なお、1Ωとは隣接両点間の線分の抵抗値。

ネット上を調べたら、(0,0)と(1,0)間の抵抗値は0.5Ωとなる解説は多いが、ほとんどが事後解釈の領域を出ない。(0,0)と(1,1)間の抵抗値が分かるのであれば、はじめて本物と認めてあげたい気持ちになる。

2. 計算法
自分の知識では全く歯が立たないので、一生懸命ネットを調べた。すごいことにそれを計算した論文をついにみつけた。

(PDF) Infinite resistive lattices

細かい中身はご自身でみてもらうことにして、必要な計算式と計算結果の一部を紹介しておく。

p11271-1.jpg
p11271-2.jpg

変数nとpはそれぞれy座標、x座標と考えてよい。対称性があるから、入れ替えても計算結果は同じになるはず。

p11271-3.jpg

上のテーブルには (5,5)までの抵抗値がリストアップされている。例えば、(0,0)-(1,1)間の抵抗値は2/Phi=0.637Ωになる。円周率Phiがここに出てくる。

ということで、論文を見る限り、数値積分をしないといけなさそう。実際にプログラムを組んで計算してみたら、上記のテーブルに近い値が出てきた。ただ、積分部分のプログラム(シンプソン法)は精度が良くないせいか、xとyの対称性に若干の誤差があったりする。UVaの自動判定システムにもAcceptedされていない。

単純そうに見えた計算問題に大変複雑な定積分が出てくるとはびっくる。漸化式では無理かな。

でも、問題文のヒントがとても気になる。もっとうまい計算法はあるような言い方。
Hint: There is a surprising Dynamic Programming solution, but how do you get it to fit under the memory requirement? 🙂 .

最後に、プログラム内に使っている、計算式に関連するコードをピックアップしておく。

int x, y;  /* 計算したい (x,y) 座標 */
double alpha(double beta)
{
    double t = cos(beta);
    return log(2-t+sqrt(3+t*(t-4)));
}
double func(double beta)
{
    double a = alpha(beta);
    val = (1.0-1/exp(x*a)*cos(y*beta))/sinh(a);
    return val;
}

上のやり方では問題ないように思われるが、確信はまだない。Acceptedされていないし。ちなみに、現時点では問題を解いたひとは二人だけ、出題者本人とUVaシステムの関係者?

新しいUVaサイトがよく落ちるね。広告を入れるのはしょうがないが、もう少しセンスがあっていいじゃないかな。無駄なスペースを取りすぎ、とくに垂直方向。肝心な問題文をスクロールしてみるという有様。IEもfirefoxも開けない。なんとかしてほしいよ。

今日から使えたGoogle CromeはSafariと同様、ちゃんとアクセスできてよかった。

#11281番の問題はよく勘違いされるようで、自分も最初では間違っていた。やはりちゃんと図形を描いてみることだね。

1. 問題文
丸い円をいくつも与えられて、さらに様々な形の三角形が与えられている。円に入る三角形を教えろ、というのが問題の意味。実際の問題では、円の代わりに孔を使ったり、三角形の代わりに立体の三角プリズムを使っているが。

2. 考え方
三角形の外接円の半径を求め、円の半径と比較すれば簡単に問題が解けそうに思われるが、大きな落とし穴が待ち構えている。

たとえば、下の三角形をよくみてみよう。ブルーの円は外接円だけど、赤の円よりは明らかに半径が大きい。つまり、鈍角三角形に対しては、その最長辺が最小円の直径になる。外接円では全然大きい。鋭角三角形に対しては、外接円が最小円になるけど。

p11281.jpg

上のことを理解すれば、問題は簡単に解けるはず。

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

つまり、自然数Nを3つの自然数 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通りしか選べないので、分割数は4になる。

p11480.png

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

  分割数 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個に分割する時の数式ってあるのだろうか。