ちなみに、以下のコードで計測すると、わたしの環境では、
$ g++ -O0 gcd.c
$ ./a.exe 1 5 10000000
seed=1, n=5, N=10000000
gcd0: 35.72 sec / 5*10000000 times
gcd1: 34.3 sec / 5*10000000 times
gcd2: 57.67 sec / 5*10000000 times
gcd3: 37.39 sec / 5*10000000 times
$ g++ -O4 gcd.c
$ ./a.exe 1 5 10000000
seed=1, n=5, N=10000000
gcd0: 28.05 sec / 5*10000000 times
gcd1: 28.52 sec / 5*10000000 times
gcd2: 29.2 sec / 5*10000000 times
gcd3: 26.28 sec / 5*10000000 times
と、なりました。前処理あるなしや、入力によるでしょうが、わたしのマシンは、割り算が遅くて、剰余演算の実装も貧弱なんでしょうね、たぶん^^
Windows XP Home Edition Ver.2002 SP2
Dell Dimension DE051 Intel(R) Pentium(R) 4 CPU 2.80GHz 1.00GB RAM
cygwin ver. 1.0.2-1
g++ v.4.1.1
お邪魔しました^^;
===============
/*
最大公約数を求める関数の 4 つの実装
0. 互除法(剰余演算子を使わない方法)
1. 互除法(剰余演算子を使う方法)
2. 互減法(普通の減算を用いる方法)…もともとの Euclid の方法
3. 互減法( 2 進表示を活用する特殊な方法)
テンプレート引数 I は、符号付き整数型
*/
template<typename I> I pp_gcd_(I *, I *);
template<typename I> I gcd0_(I, I);
template<typename I> I gcd1_(I, I);
template<typename I> I gcd2_(I, I);
template<typename I> I gcd3_(I, I);
/* 互除法(剰余演算子を使わない方法) */
template<typename I> I gcd0(I a, I b)
{
I c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd0_(a, b);
else c = gcd0_(b, a);
}
return c;
}
/* 互除法(剰余演算子を使う方法) */
template<typename I> I gcd1(I a, I b)
{
I c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd1_(a, b);
else c = gcd1_(b, a);
}
return c;
}
/* 互減法(通常の減法を使う方法) */
template<typename I> I gcd2(I a, I b)
{
I c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd2_(a, b);
else c = gcd2_(b, a);
}
return c;
}
/* 互減法(2進表現を使う方法) */
template<typename I> I gcd3(I a, I b)
{
I c = pp_gcd_(&a, &b);
if (c < 0) c = gcd3_(a, b);
return c;
}
/*
前処理:
・双方 0 なら 0 を返す
・一方が 0 ならもう一方の絶対値を返す
・絶対値に設定しなおし、双方非 0 で、絶対値が同じなら、その絶対値を返す
・絶対値に設定しなおし、双方非 0 で、絶対値が違うなら -1 を返す
*/
template<typename I> I pp_gcd_(I *a, I *b)
{
if (!*a) {
return *b < 0 ? -*b : *b;
} else if (!*b) {
return *a < 0 ? -*a : *a;
}
if (*a < 0) *a = -*a;
if (*b < 0) *b = -*b;
if (*a == *b) return *a;
return -1;
}
/* a > b > 0 を仮定(a > 0 && b > 0 で十分であるが) */
template<typename I> I gcd0_(I a, I b)
{
while (1) {
if (!(a -= a / b * b)) return b;
if (!(b -= b / a * a)) return a;
}
}
/* a > b > 0 を仮定(a > 0 && b > 0 で十分であるが) */
template<typename I> I gcd1_(I a, I b)
{
while (1) {
if (!(a %= b)) return b;
if (!(b %= a)) return a;
}
}
/* a > b > 0 を仮定 */
template<typename I> I gcd2_(I a, I b)
{
while (1) {
do { a -= b; } while (a > b);
if (a == b) return b;
do { b -= a; } while (b > a);
if (a == b) return a;
}
}
/* a > 0 && b > 0 を仮定 */
template<typename I> I gcd3_(I a, I b)
{
/* 2の何乗で割り切れるか見ておく */
I p2a = 0, p2b = 0;
while (!(a & 1)) { a >>= 1; ++p2a; }
while (!(b & 1)) { b >>= 1; ++p2b; }
/* a も b も奇数である。a == b になるまで引き算する(a と b の差は偶数) */
while (a != b) {
if (a < b) {
b -= a;
while (!(b & 1)) b >>= 1;
} else {
a -= b;
while (!(a & 1)) a >>= 1;
}
}
return a << (p2a < p2b ? p2a : p2b); /* 2のべき乗因子をかけて返す */
}
template<typename I> void measure(unsigned long, I, I);
#include <iostream>
#include <sstream>
// 前処理が必要な場合、定義する
/* #define NEED_PP_GCD_ */
// 非標準の random() を使う場合、定義する
#define USE_RANDOM_
int main(int argc, char *argv[])
{
typedef int int_type;
unsigned long seed = 1;
int_type n = 0, N = 0;
if (argc > 1) { std::istringstream iss(argv[1]); iss >> seed; }
if (argc > 2) { std::istringstream iss(argv[2]); iss >> n; }
if (argc > 3) { std::istringstream iss(argv[3]); iss >> N; }
if (n <= 0) n = 100;
if (N <= 0) N = 1000000;
std::cerr << "seed=" << seed << ", n=" << n << ", N=" << N << '\n';
measure(seed, n, N);
}
#include <ctime>
/*
ストップウォッチ(CPU 時間)
*/
class StopClock
{
std::clock_t total; // CPU 時間合計
std::clock_t base; // 計測開始時刻
bool running; // 計測中か否か
public:
// 引数が true なら、構築後すぐに計測状態になる。
StopClock(bool start_now = false)
{
reset();
if (start_now) start();
}
// 計測開始。ただし、すでに計測中なら何もしない。
void start()
{
if (running) return;
running = true;
base = std::clock();
}
// 計測を終了し、合計する。ただし、計測外なら何もしないが、できるだけ精確に
// 時間を計測したいため、計測中か否かを調べずにすぐ時刻を取得している。
void stop()
{
std::clock_t t;
t = std::clock();
if (!running) return;
total += t - base;
running = false;
}
// 合計時間をゼロにして、計測外の状態にする。
void reset()
{
total = 0;
running = false;
}
// 合計時間を返す。ただし、計測中なら計測を終了させ、合計時間を返す。
double time()
{
stop();
return total / static_cast<double>(CLOCKS_PER_SEC);
}
};
#include <sys/time.h>
/*
ストップウォッチ(経過時間)
*/
class StopWatch
{
timeval total; // 経過時間合計
timeval base; // 計測開始時刻
bool running; // 計測中か否か
static const suseconds_t one_second = 1000000; // 一秒は百万マイクロ秒
public:
// 引数が true なら、構築後すぐに計測状態になる。
StopWatch(bool start_now = false)
{
reset();
if (start_now) start();
}
// 計測開始。ただし、すでに計測中なら何もしない。
void start()
{
if (running) return;
running = true;
gettimeofday(&base, 0);
}
// 計測を終了し、合計する。ただし、計測外なら何もしないが、できるだけ精確に
// 時間を計測したいため、計測中か否かを調べずにすぐ時刻を取得している。
void stop()
{
timeval t;
gettimeofday(&t, 0);
if (!running) return;
time_t s = t.tv_sec - base.tv_sec;
suseconds_t m = t.tv_usec;
if (t.tv_usec >= t.tv_usec) {
m -= base.tv_usec;
} else {
--s;
m += one_second - base.tv_usec;
}
total.tv_sec += s;
total.tv_usec += m;
if (total.tv_usec >= one_second) {
total.tv_usec -= one_second;
++total.tv_sec;
}
running = false;
}
// 合計時間をゼロにして、計測外の状態にする。
void reset()
{
total.tv_sec = total.tv_usec = 0;
running = false;
}
// 合計時間を返す。ただし、計測中なら計測を終了させ、合計時間を返す。
double time()
{
stop();
return total.tv_sec + static_cast<double>(total.tv_usec) / one_second;
}
};
#ifdef USE_RANDOM_
#include <stdlib.h>
#else
#include <cstdlib>
#endif
template<typename I> void f0(I N, const I a[], const I b[], I c[])
{
#ifdef NEED_PP_GCD_
for (I i = 0; i < N; ++i) c[i] = gcd0(a[i], b[i]);
#else
for (I i = 0; i < N; ++i) c[i] = gcd0_(a[i], b[i]);
#endif
}
template<typename I> void f1(I N, const I a[], const I b[], I c[])
{
#ifdef NEED_PP_GCD_
for (I i = 0; i < N; ++i) c[i] = gcd1(a[i], b[i]);
#else
for (I i = 0; i < N; ++i) c[i] = gcd1_(a[i], b[i]);
#endif
}
template<typename I> void f2(I N, const I a[], const I b[], I c[])
{
#ifdef NEED_PP_GCD_
for (I i = 0; i < N; ++i) c[i] = gcd2(a[i], b[i]);
#else
for (I i = 0; i < N; ++i)
if (a[i] > b[i]) c[i] = gcd2_(a[i], b[i]);
else if (a[i] < b[i]) c[i] = gcd2_(b[i], a[i]);
else c[i] = a[i];
#endif
}
template<typename I> void f3(I N, const I a[], const I b[], I c[])
{
#ifdef NEED_PP_GCD_
for (I i = 0; i < N; ++i) c[i] = gcd3(a[i], b[i]);
#else
for (I i = 0; i < N; ++i) c[i] = gcd3_(a[i], b[i]);
#endif
}
#include <cmath>
#include <iomanip>
// iostream フォーマット番兵
template< typename Ch, typename Tr = std::char_traits<Ch> >
struct BasicIosFmtSentry {
std::basic_ios<Ch, Tr> *pbi;
std::basic_ios<Ch, Tr> store;
BasicIosFmtSentry(std::basic_ios<Ch, Tr> *pbi_) : pbi(pbi_), store(0)
{ store.copyfmt(*pbi); }
~BasicIosFmtSentry() { pbi->copyfmt(store); }
};
typedef BasicIosFmtSentry<char> IosFmtSentry;
template<typename I> void measure(unsigned long seed, I n, I N)
{
I *mem = new I[3 * N * sizeof(I)];
I *a = mem, *b = mem + N, *c = mem + 2 * N;
#ifdef USE_RANDOM_
srandom(seed);
#else
std::srand(static_cast<unsigned int>(seed));
#endif
for (I i = 0; i < N; ++i) {
#ifdef USE_RANDOM_
a[i] = static_cast<I>(random());
b[i] = static_cast<I>(random());
#else
a[i] = static_cast<I>(rand());
b[i] = static_cast<I>(rand());
#endif
#ifndef NEED_PP_GCD_
if (!a[i]) a[i] = 1;
if (!b[i]) b[i] = 1;
#endif
}
std::streamsize prec = static_cast<std::streamsize>(
1 + std::log10(static_cast<double>(CLOCKS_PER_SEC))
);
IosFmtSentry ifs(&std::cout);
std::cout << std::setprecision(prec) << std::showbase;
StopClock sc;
sc.reset(); sc.start();
for (I i = 0; i < n; ++i) f0(N, a, b, c);
sc.stop();
std::cout << "gcd0: " <<
sc.time() << " sec / " << n << '*' << N << " times \n";
sc.reset(); sc.start();
for (I i = 0; i < n; ++i) f1(N, a, b, c);
sc.stop();
std::cout << "gcd1: " <<
sc.time() << " sec / " << n << '*' << N << " times\n";
sc.reset(); sc.start();
for (I i = 0; i < n; ++i) f2(N, a, b, c);
sc.stop();
std::cout << "gcd2: " <<
sc.time() << " sec / " << n << '*' << N << " times\n";
sc.reset(); sc.start();
for (I i = 0; i < n; ++i) f3(N, a, b, c);
sc.stop();
std::cout << "gcd3: " <<
sc.time() << " sec / " << n << '*' << N << " times\n";
delete[] mem;
}