自然乗数演算に限られ,ビット演算法則を用いてC++に組み込まれたpow()関数を迅速に計算する方法を理解する.

TL;DR

double pow_int(double base, int exp)
{
    double result = 1.;
    while (exp)
    {
        if (exp & 1)
            result *= base;
        exp >>= 1;
        base *= base;
    }
    return result;
}

n/a原理

実数aaaと自然数nnnについて,ana^nanの計算を見てみよう.
ここで、nnnは、バイナリ変換によって次のように表すことができる.
n=∑i=0bi2i, where bi=0 or 1n =\sum_{i=0} b_i 2^i\text{, where }b_i=0\text{ or }1n=∑i=0​bi​2i, where bi​=0 or 1
上記の式を使用すると、ana^nanは次の変更を行います.
\begin{equation}\begin{aligned} a^n & = a^{\sum_{i=0} b_i 2^i}\\\\& =\prod_{i=0} {a^{b_i 2^i}}\text{ (}\because x^{m+n} = x^m x^n\text{)}\\\\& =\prod_{i=0} ({a^{2^i}})^{b_i}\text{ (}\because x^{mn} = (x^m)^n\text{)}\\\\& = (b_0\times a)\times (b_1\times a^2)\times (b_2\times (a^2)^2)\times (b_3\times ((a^2)^2)^2)\times ...\end{aligned}\end{equation}
この式の意味を説明します.

nnをバイナリで表す場合、低いビット数から決定する(b 0、b 1、...b 0、b 1、...b 0、b 1、...)

各bib ibi,a,a 2,(a 2)2,...a, a^2, (a^2)^2, ...a,a2,(a2)2,...の形式で平方を続けます

bib ibiが0ならそのままスキップして、
bib ibiが1の場合、2から得られたa 2 ia^{2^i}a 2 iを加算する.

例

353^535を計算します.ここで、a=3 a=3 a=3、n=5 n=5 n=5である.

まずnnnをバイナリに変換してbib ibiを求める.
n=5=1+4=1×20+0×21+1×22n=5=1+4=1\times2^0+0\times2^1+1\times2^2n=5=1+4=1×20+0×21+1×22
i.e. b0=1, b1=0, b2=1, and bk=0 for k>=3\text{i.e. }b_0=1\text{, } b_1=0\text{, } b_2 = 1\text{, and } b_k = 0\text{ for }k>=3i.e. b0​=1, b1​=0, b2​=1, and bk​=0 for k>=3

次に、低い桁数からチェックを開始し、乗算を実行します.積算積の初期値は、積の一定因子1に設定されます.

aaaは3です.
b 0 b 0 b 0 b 0を確認します.b 0 b 0 b 0 b 0は111なので、1に3を乗じる.
結果は3です.

a 2 a^2 a 2は9です.
b 1 b 1 b 1 b 1を確認します.b 1 b 1 b 1 b 1は000で、何も乗じません.
結果は依然として3です.

(a 2)2(a^2)^2(a 2)2は81です.
b2b2 2b2.b 2 b 2 b 2 b 2は111であり、3に81を乗じる.
結果は243であった.

以降のbkb kbkはいずれも000であるため,最終結果は243であった.

コードに適用(C++)

上記の内容をコードに変換します.
変数名の意味を考慮すると、上記の式ではaaaはbase、nnnはexpと命名される.

double pow_int(double base, int exp)
{
    double result = 1.;
    while (exp)
    {
        if (exp & 1)
            result *= base;
        exp >>= 1;
        base *= base;
    }
    return result;
}

result変数の初期値を乗算の定式化1に設定し、積算乗算を実行します.

値がexpの場合、while繰返しが実行されます.

expをバイナリで表す場合の最低桁数を求めます.( exp & 1 )
この値が1の場合、bib ibiは1です.したがって、現在の値baseはresultに加算される.
この値が0の場合、bib ibiは0です.だからresultに何も乗じない
例えば、expのb 0 b 0 b 0 b 0 b 0値が1である場合、このif条件は、最初の繰り返し文の実行時に合致し、result *= base;構文が、上記の説明のaaaに乗算されるのと同じ動作を実行する.

expから右に1つのビット切り替えを実行します.( exp >>= 1; )
これにより、さっき表示したbib ibiを破棄し、対応するbi+1 b{i+1}bi+1を最低位置に移動するプロセスを繰り返します.

baseを掛けます.( base *= base; )
これによって(…(a 2)2)…)2(...(a^2)^2)...)^2(...(a2)2)...)2に対応する繰返し演算が実行されます.

最終的に累加乗算を実行した結果を返します.( return result; )

pow()関数と比較するためのサンプルコード

コード#コード#

#include <iostream>

double pow_int(double base, int exp)
{
    double result = 1.;
    while (exp)
    {
        if (exp & 1)
            result *= base;
        exp >>= 1;
        base *= base;
    }
    return result;
}

int main()
{
    const int N = 100000000;
    double base = 12.345;
    int exp = 67;
    double result1, result2;

    std::cout << base << "^" << exp << "\n\n";

    auto startT1 = clock();
    for (int i = 0; i < N; ++i)
        result1 = pow(base, exp);
    auto endT1 = clock();
    std::cout << "Calculated by pow() function\n";
    std::cout << "Answer: " << result1 << "\n";
    std::cout << "Elapsed time: " << (endT1 - startT1) << " ms\n\n";
    
    auto startT2 = clock();
    for (int i = 0; i < N; ++i)
        result2 = pow_int(base, exp);
    auto endT2 = clock();
    std::cout << "Calculated by pow_int() function\n";
    std::cout << "Answer: " << result2 << "\n";
    std::cout << "Elapsed time: " << (endT2 - startT2) << " ms\n\n";

    return 0;
}

結果

3.00 GHzベースのi 5-9500