上式の通り、ごく簡単な、ベキ級数の合計を算出する問題。変数 l, h, k の値はともに1～15000000、hとlとの差は1000以内。

kの値が10以内なら、探せば、それなりの計算式が出てくるが、1500万ではそんな式は存在しないだろう。

力任せで、C言語の数学関数をそのまま使うと、無限大になってしまい、計算不能になる。

考え方
　合計値は級数の後ろの項によってほぼ決まるので、級数の末尾から計算することにしよう。

大事なことは、個々の項の値を、仮数部と指数部に分けておく。

例えば、999¹⁰⁰、998¹⁰⁰の値を
　　999¹⁰⁰　仮数部 0.9047921471　指数部 e+300
　　998¹⁰⁰　仮数部 0.8185668046　指数部 e+300
に分ける。

指数部の値が1つ違うと、項の値として1桁も違ってくる。指数部の値が5も違うと、5桁も右の部分にしか影響しないので、その後の計算を打ち切っても全体の合計値にほとんど影響がないだろうと考えてよい。例えば、
　　880¹⁰⁰　仮数部 0.2807160311　指数部 e+295
当たりから、その後の合計計算を打ち切る。基数やkの値が大きければ大きいほどすぐに打ち切ることになる。

また、合計の計算では、級数の末尾から項ひとつひとつ足していくが、個々の項の指数部は、末尾項の指数部の値に（つまり、最大指数部）に合わせて、仮数部の小数点をシフトさせる。
　上の例では、級数の末尾項が999¹⁰⁰なら、個々の項の指数部の値をe+300に統一して、仮数部の値を合計する。
　　970¹⁰⁰　仮数部 0.0475525075　指数部 e+300

なお、シフトする前の仮数部、指数部は以下の数学関数でも求めることができる。

　ベキ乗 x^k → 仮数部 pow(10, modf(k*log10(x), ＆e))　指数部 左の e

最後に、元問題の検証データを載せておく（入力データの各行は左から、l, h, kの値）

10 10 1000000
14999001 15000100 1
14999001 15000000 15000000
0 0 15000000
0 1 15000000
0 0 1

Case 0001: 0.100000e0001000001
Case 0002: 0.164995e0000000011
Case 0003: 0.121628e0107641370
Case 0004: 0.000000e0000000001
Case 0005: 0.100000e0000000001
Case 0006: 0.000000e0000000001

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