ユークリッドの除去法アルゴリズム

最大公約数を求める際ユークリッドの除去法を使ったアルゴリズムを考える場合、計算量はO(log max{x,y})となる理由を教えて下さい。

簡単な擬似コードも教えてもらえるとありがたいです。

No.1 回答者： mikaemi 回答日時： 2007/09/03 23:03

数学の問題かとも思いますが^^

「ラメの定理」として知られていますね。

あるいは、もう少し精密にして、「φ=(√5+1)/2 とすると（フィボナッチ数列で出てくる数値）、互除法は、φを底にした対数関数log(min(m,n))以上の最小の整数回の除法で完了する」となるようです。

たとえば、一松信 著『代数学入門第一課』近代科学社 に載ってます。

ラメの定理（Lame's theorem）で検索したら見つかると思いますよ、たぶん^^ 英語かもしれないですが。

No.2 回答者： asuncion 回答日時： 2007/09/03 23:17

もう一方の問題を含めて、何となく

「情報処理基礎」あたりの前期試験問題っぽい…。

No.3 回答者： mikaemi 回答日時： 2007/09/03 23:18

証明はそんなに難しくも無く、長くも無いので、証明を書こうとしたんですが、数式を生のテキストでうまく書くのが面倒なのでやめてしまいました^^; 生のテキストで書いても読みづらいですし^^

一応、http://www.cs.fsu.edu/~cop4531/slideshow/chapter …　にそれっぽい証明があるようです。

No.4 回答者： mikaemi 回答日時： 2007/09/03 23:19

あぁ、なるほど。

試験問題とか課題の質問か（笑）

No.5 回答者： Oh-Orange 回答日時： 2007/09/04 02:14

★回答ではないが気になったので１つ。

・ユークリッドの除去法(じょきょほう)⇒聞いた事がないよ。
　ユークリッドの互除法(ごじょほう)⇒聞いた事がある。分数計算でお世話になった。

参考文献:
　http://ja.wikipedia.org/wiki/%E3%83%A6%E3%83%BC% …
　http://www.rimath.saitama-u.ac.jp/lab.jp/fsakai/ …
　http://www2.cc.niigata-u.ac.jp/~takeuchi/tbasic/ …
・以上。

No.6 回答者： mikaemi 回答日時： 2007/09/04 06:20

まったく、回答ではないですが（笑）ちなみに、Leftymouse さんのマシン環境で、以下の４つ（どれも結局、余りを求めているが^^）のうちのどれが一番高速ですか？

＝＝＝＝
/*
最大公約数を求める関数の 4 つの実装

0. 互除法(剰余演算子を使わない方法)
1. 互除法(剰余演算子を使う方法)
2. 互減法(普通の減算を用いる方法)…もともとの Euclid の方法
3. 互減法( 2 進表示を活用する特殊な方法)

int_type は、符号付き整数型
*/

typedef int int_type;
/* typedef long long int_type; */

int_type pp_gcd_(int_type *, int_type *);
int_type gcd0_(int_type, int_type);
int_type gcd1_(int_type, int_type);
int_type gcd2_(int_type, int_type);
int_type gcd3_(int_type, int_type);

/* 互除法(剰余演算子を使わない方法) */
int_type gcd0(int_type a, int_type b)
{
int_type c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd0_(a, b);
else c = gcd0_(b, a);
}
return c;
}

/* 互除法(剰余演算子を使う方法) */
int_type gcd1(int_type a, int_type b)
{
int_type c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd1_(a, b);
else c = gcd1_(b, a);
}
return c;
}

/* 互減法(通常の減法を使う方法) */
int_type gcd2(int_type a, int_type b)
{
int_type c = pp_gcd_(&a, &b);
if (c < 0) {
if (a > b) c = gcd2_(a, b);
else c = gcd2_(b, a);
}
return c;
}

/* 互減法(2進表現を使う方法) */
int_type gcd3(int_type a, int_type b)
{
int_type c = pp_gcd_(&a, &b);
if (c < 0) c = gcd3_(a, b);
return c;
}

/*
前処理:
・双方 0 なら 0 を返す
・一方が 0 ならもう一方の絶対値を返す
・絶対値に設定しなおし、双方非 0 で、絶対値が同じなら、その絶対値を返す
・絶対値に設定しなおし、双方非 0 で、絶対値が違うなら -1 を返す
*/
int_type pp_gcd_(int_type *a, int_type *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 で十分であるが) */
int_type gcd0_(int_type a, int_type b)
{
while (1) {
if (!(a -= a / b * b)) return b;
if (!(b -= b / a * a)) return a;
}
}

/* a > b > 0 を仮定(a > 0 && b > 0 で十分であるが) */
int_type gcd1_(int_type a, int_type b)
{
while (1) {
if (!(a %= b)) return b;
if (!(b %= a)) return a;
}
}

/* a > b > 0 を仮定 */
int_type gcd2_(int_type a, int_type 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 を仮定 */
int_type gcd3_(int_type a, int_type b)
{
/* 2の何乗で割り切れるか見ておく */
int_type 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のべき乗因子をかけて返す */
}

#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[])
{
if (argc < 3) {
fprintf(stderr, "Usage: %s integer integer\n", argv[0]);
return 1;
}
int a = atoi(argv[1]), b = atoi(argv[2]);
printf("gcd0(%d, %d)=%d\t", a, b, gcd0(a, b));
printf("gcd1(%d, %d)=%d\t", a, b, gcd1(a, b));
printf("gcd2(%d, %d)=%d\t", a, b, gcd2(a, b));
printf("gcd3(%d, %d)=%d\n", a, b, gcd3(a, b));
/*
long long a = atoll(argv[1]), b = atoll(argv[2]);
printf("gcd0(%lld, %lld)=%lld\t", a, b, gcd0(a, b));
printf("gcd1(%lld, %lld)=%lld\t", a, b, gcd1(a, b));
printf("gcd2(%lld, %lld)=%lld\t", a, b, gcd2(a, b));
printf("gcd3(%lld, %lld)=%lld\n", a, b, gcd3(a, b));
*/
return 0;
}

No.7 ベストアンサー 回答者： mikaemi 回答日時： 2007/09/04 10:50

ちなみに、以下のコードで計測すると、わたしの環境では、

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

No.8 回答者： mikaemi 回答日時： 2007/09/04 11:09

あっ。

。。main() 中の typedef は、

　typedef long int_type;

じゃないと sizeof(int) < sizeof(long) のマシンだとダメですね。失礼しました^^

No.9 回答者： mikaemi 回答日時： 2007/09/04 11:27

あと、

　std::cout << std::setprecision(prec) << std::showbase;
じゃなくて、
　std::cout << std::setprecision(prec) << std::showpoint;
でした^^;

No.10 回答者： mikaemi 回答日時： 2007/09/04 11:41

最後に、型変更によるオーバフローを考慮して、

　#ifndef NEED_PP_GCD_
　if (!a[i]) a[i] = 1;
　if (!b[i]) b[i] = 1;
　#endif

じゃなくて、

　#ifndef NEED_PP_GCD_
　if (a[i] <= 0) a[i] = 1;
　if (b[i] <= 0) b[i] = 1;
　#endif

のようにしておくべきでしたね^^ ではでは。