プログラミングコンテストの数学

Topcoder SRM等の問題を見ると，ときおり数学的知識を要する問題が出てきます。
学業から離れて数年も経つと，仕事に関係しない知識はどんどんこぼれるもので，不意に数学の知識を問う問題に出会うと戸惑います。

この辺りのところを，今後のコンテストのためにメモにまとめておこうと思います。

素数

まずは素数判定から。素数判定には様々な手法が提案されていますが，小さな数の素数判定であれば次の単純な判定処理で十分ではないかと思います。

bool isPrime(unsigned int n)
{
    if( n == 2 ) return true;
    if( (n&1) == 0 || n == 1 ) return false;// 1と，2以外の偶数は素数でない
    int upper = (int)sqrt((double)n);
    for( int i = 3; i <= upper; i += 2 ){// 3以上の奇数を評価
        if( n % i == 0 ) return false;
    }
    return true;
}

判定対象nより小さい数で割り切れるかを見ます。割り切れた場合は素数でありません。割り切れるかの判定は&radicnまでで十分です。
ループ回数が入力値の平方根/2ですので，UINT_MAX(=4G)であってもたかだか30000回ほどのループ回数で済みます。

---

素数のリストが欲しい場合は，「エラストテネスのふるい」を用いるのが一般的です。

void sieve( int n, vector<int>& primes )
{
    if( n < 2 ) return;
    primes.push_back(2);
    if( n == 2 ) return;

    vector<bool> p(n/2,true);// p[x] => 2*x+3が素数であるかを示す
    int upper = (int)sqrt((double)n);
    upper = (upper-3)/2;

    for( int i = 0; i <= upper; i++ ){
        int t = i*2+3;
        if( p[i] ){
            primes.push_back(t);
            p[i] = false;
            for( int j = (t+1)*i; j < n/2; j += t ){
                p[j] = false;
            }
        }
    }
    for( int i = upper; i < (n-1)/2; i++ )
        if( p[i] )
            primes.push_back(i*2+3);
}

このコード例では，nまでの素数をvector primes に格納します。

素数判定用のワーク配列 p を n/2 のサイズで用意し，素数であることを示すtrueで初期化します。コード例では，配列 p を3以上の奇数を示すインデックス値の配列として考えています。
「エラストテネスのふるい」では次の処理を繰り返します。

1. 配列 p の中で，trueである最も小さい数を primes へ入れる
2. 項1で選択した数の倍数をfalseにする

ふるいをかける上限は，素数の判定と同様に &radicn までで十分です。また，素数 p_i の倍数をfalseにする部分では p_i×p_i 以降の数について処理すれば十分となります。これ以前の数は p_i より小さい素数の倍数処理で既に false にされています。

最大公約数

規約分数を求めるときに最大公約数が必要となります。 最大公約数は，ユークリッドの互助法で求められます。

unsigned int gcd( unsigned int a, unsigned int b )
{
    unsigned int t;
    if( a < b ){
        t = a; a = b; b = t;
    }
    while( b != 0 ){
        t = a % b;
        a = b;
        b = t;
    }
    return a;
}

ユークリッド互助法は，下記の式が( b≠0 のとき)成り立つことを使っています。

この関係式を繰り返し適用することで最小公倍数を取得します。

３個以上の数の最大公約数では次の関係式が成り立ちます。

---

参考情報として「最小公倍数」は次の式で求められます。

unsigned int lcm( unsigned int a, unsigned int b )
{
    return a*b/gcd(a,b);
}

基数変換

10進数と n進数の相互変換を考えます。10進数 dと n進数値の次の関係式を用います。

▼n進数への変換

void toNAdic( unsigned int decimal, vector<int>& coef, int base )
{
    do{
        unsigned int c = decimal % base;
        coef.push_back(c);
        decimal /= base;
    }while( decimal != 0 );
}

n進数の各桁の数値がvector coef に格納されます。

▼10進数への変換

unsigned int toDecimal( const vector<int>& coef, int base )
{
    unsigned int result = 0, u = 1;
    vector<int>::const_iterator it;
    for( it = coef.begin(); it != coef.end(); it++ ){
        result += *it * u;
        u *= base;
    }
    return result;
}

n進数の各桁の数値vector coef を与えると，10進の数値が返ります。

組み合わせ・順列

憶えていそうであやふやになりがちな「組み合わせ」と「順列」を整理しておきたいと思います。

組み合わせ
順列

▼組み合わせ

unsigned int combination( unsigned int n, unsigned int r )
{
    unsigned int numerator = 1, denominator = 1;
    if( n-r < r ) r = n-r;
    if( r == 0 )  return 1;
    for( unsigned int i = n; i > (n-r); i-- )
        numerator *= i;
    for( unsigned int i = 2; i <= r; i++ )
        denominator *= i;
    return numerator/denominator;
}

[long long版]

unsigned long long combination( unsigned long long n, unsigned long long r )
{
    unsigned long long p = 1, a= 1;
    if( n-r < r ) r = n-r;
    if( r == 0 )  return 1;
    for( unsigned long long i = n; i > (n-r); i-- )
        p *= i;
    for( unsigned long long i = 2; i <= r; i++ )
        a *= i;
    return p/a;
}

▼順列

unsigned int permutation( unsigned int n, unsigned int r )
{
    unsigned int result = 1;
    for( unsigned int i = n; i > (n-r); i-- )
        result *= i;
    return result;
}

[long long版]

unsigned long long permutation( unsigned long long n, unsigned long long r )
{
    unsigned long long result = 1;
    for( unsigned long long i = n; i > (n-r); i-- )
        result *= i;
    return result;
}

どちらのコード例もオーバーフローを考慮していませんので，適切な範囲で値を渡してください。

参考までに，階乗計算でオーバフローせずに計算可能な値は次の値までとなります。

long :	12! = 479,001,600
long long :	20! = 2,432,902,008,176,640,000

三角関数

C/C++における三角関数の仕様を纏めておきます。

三角関数	引数	戻り値
sin(r)	ラジアン	sin値: -1 <= ret <= 1
cos(r)	ラジアン	cos値: -1 <= ret <= 1
tan(r)	ラジアン	tan値: π/2で∞であるが，実際には近似値を使うため値が返る。
asin(t)	-1 <= t <= 1	ラジアン: -π/2 <= ret <= π/2
acos(t)	-1 <= t <= 1	ラジアン: 0 <= ret <= π
atan(t)	任意	ラジアン: -π/2 <= ret <= π/2
atan2(y, x)	任意	ラジアン: -π <= ret <= π　　[x=0かつy=0のときは，0が返る]

モジュロ演算

モジュロ演算（剰余演算）では，次の式が成立します。

足し合わせた結果の剰余演算は，それぞれの要素の剰余演算の結果を足し合わせたものと同じ値となります。

グレイ符号・ビットカウント

最後はあまり数学的なトピックではありませんが，グレイ符号とビットカウントについて記しておきます。過去のSRMで調べておかなければと思った項目です。

▼グレイ符号
グレイ符号はビット符号化方式であり，隣り合う値がハミング距離で１しか離れない符合となります。詳細はWikipediaに譲ります。

unsigned int GrayCode(unsigned int val)
{
    return val^(val>>1);
}

▼ビットカウント
数値の2進表現において，1のビットがいくつあるかを数え上げます。cell(spu)だと，cntbですね。

検索してみると次のコードが見つかります。

unsigned int cntb(unsigned int val)
{
    val =  (val & 0x55555555) + (val >> 1 & 0x55555555);
    val =  (val & 0x33333333) + (val >> 2 & 0x33333333);
    val =  (val & 0x0f0f0f0f) + (val >> 4 & 0x0f0f0f0f);
    val =  (val & 0x00ff00ff) + (val >> 8 & 0x00ff00ff);
    return (val & 0x0000ffff) + (val >>16 & 0x0000ffff);
}

最初は1ビット同士の足し算，その結果が2ビットで表されるので続いて2ビット同士の足し算を行っています。これを繰り返し，最後は16ビットの足し算を行って1のビットの合計が求まります。

参考文献

・C言語による最新アルゴリズム事典
・Mathematics for TopCoders : Topcoder Algorithm Tutorial