binary exponentiation

任意の数の累乗を計算する時には、for文で累乗する回数だけ数を掛け合わせる、というのは小学生でも思いつく。だが、累乗する回数が大きくなったとき、例えば1億乗とかやられたときには極端に遅くなるのは明らかだ。

二進累乗法というやり方がある。これは2の累乗に乗数を展開することで、掛ける回数をビット数＋数回に抑えるという考え方だ。累乗する数を1ビットシフトして掛けることで、x^2, x^4, x^8, x^16 ... といった具合で累乗する回数を2の累乗倍に増やしていくことで、累乗する回数を減らすことができる。

累乗する回数を減らすことは、累乗する対象の値が double だった場合に誤差を減らす効果をもたらす。
.... って、言葉で言ってもわけわかんないすよね(´ー｀; )

	#include <stdio.h>
	/*
	* param の num 乗を出力する関数(二進累乗法の実装)
	*
	* TODO: オーバーフロー/アンダーフローのチェックをしていない
	* TODO: num が負の時に対応していない
	*/
	static void power(double param, long long int num)
	{
	double result = 1;
	double power_base = param;
	if (param == 0) {
	printf("0\n");
	return;
	}
	if (num == 0) {
	printf("1\n");
	return;
	}
	while (num != 0) {
	//
	// 乗数の最下位ビットをチェックし、1ならば
	// 最終結果に今までの param^2x 分を掛ける
	// これは以下の役割を果たす
	//
	// - num が奇数であった場合に、result の
	// 初期状態をpower_base に等しくする
	// （つまり、一回掛けた状態にする)
	// - 最下位ビットが1の場合をとらえて、result
	// に今までためてきた power_base^2x の値
	// を掛ける。この回数が少ないほど誤差が少なく
	// なる。
	//
	if ((num & 1) == 1) {
	result *= power_base;
	}
	//
	// 乗数を1ビットシフトし、
	// power_base^2 する。これにより、以下の
	// ような形になる
	//
	// power_base^2 -> power_base^4 -> power_base^8
	// power_base^16 -> power_base^32 ...
	//
	// これにより、乗数のビット数 + 数回の掛け算により
	// 大きな乗数でも値を計算できるようになる
	//
	num = num >> 1;
	power_base *= power_base;
	}
	printf("%.20f\n", result);
	}
	int main(int argc, char *argv[])
	{
	double base = 0;
	long long int num = 0;
	// no error processing
	printf("input base number:");
	scanf("%lf", &base);
	printf("input multiplier:");
	scanf("%Ld", &num);
	power(base, num);
	return 0;
	}

view raw power.c hosted with ❤ by GitHub

[ Update September 11th 2012 06:32 JST by m ]

べき乗のアルゴリズムのうち、上位桁から計算する方法 を上記はコードにしたものだ。これの計算量は 最大でも 2 * logN (N は累乗する回数) なので、O(logN) となる。