8 способов нахождения наибольшего общего делителя
В адаптации на C++ код выглядел примерно так:
#include <iostream>
using namespace std;
int main() {
int a, b;
cin >> a >> b;
for (int i = a; i > 0; i--) {
if (a % i == 0 && b % i == 0) {
cout << "nod = " << i;
break;
}
}
return 0;
}
Подготовка
Было решено:
- перейти от типа
intк типуlong, что бы иметь больший простор для маневра; - оформить вычисление НОД в виде функции;
- в функции
mainвызывать альтернативные версии функции вычисления НОД и сравнивать их производительность.
Очевидный прототип функции: long gcd(long a, long b). Имя функции от англ. Greatest Common Divisor – т.е. НОД.
Для такой несложной задачи я не стал связываться с профилировщиком, а использовал примитивный тайминг (подробности можно посмотреть в статье «Оптимизация кода через ручной тайминг»).
Окончательный вариант «испытательного стенда» приведён в конце статьи. А в процессе исследований я пользовался его упрощённым вариантом, обеспечивавшим запуск одной из испытуемых функций и таймирование.
В коде я не пользоваться библиотечными функциями, что бы максимальный объём кода был контролируемым. Использование библиотечных функций типа min или swap — отдельная тема для экспериментов. Оставляю её для особо дотошных читателей.
01 перебор от произвольного числа
long gcd01(long a, long b) {
long nod = 1L;
for (long i = a; i > 0; i--) {
if (a % i == 0 && b % i == 0) {
nod = i;
break;
}
}
return nod;
}
Этот алгоритм – стартовая точка исследования. Идея алгоритма очень проста: гоним переменную цикла от первого числа до 1. Если оба числа делятся на переменную цикла без остатка, значит переменная цикла и равна НОД; цикл можно завершить досрочно. Если цикл прошёл до конца, значит для этих чисел НОД равен 1.
Очевидный недостаток – несимметричность относительно аргументов. Очевидно, что НОД меньше или равен меньшему из двух чисел. Поэтому гнать цикл от большего числа не имеет смысла.
02 перебор от минимального числа
long min(long a, long b) {
return a > b ? b : a;
}
long gcd02(long a, long b) {
long nod = 1L;
for (long i = min(a, b); i > 0; i--) {
if (a % i == 0 && b % i == 0) {
nod = i;
break;
}
}
return nod;
}
Просто добавляем простейшую функцию для вычисления минимального числа для пары чисел и инициализируем переменную цикла меньшим из двух чисел. В половине случаев такая оптимизация работать не будет (когда первый аргумент и так меньше второго); зато в другой половине случаев выигрыш по времени может быть весьма значительным.
03 алгоритм с разложением на делители
Но второй вариант так и остался алгоритмом с последовательным перебором. Что можно попробовать ещё?
Из школьного курса математики известно, что НОД можно найти из разложения чисел на простые множители.
a = m1 * m2 * m3 * ... * mN где m – простые числа
НОД будет равен произведению простых множителей, общих для обоих чисел. Например:
24 = (2) * 2 * 2 * (3) 54 = (2) * (3) * 3 * 3 общие множители выделены скобками НОД(24, 54) = 2 * 3 = 6
Реализуем эту идею:
long gcd03(long a, long b) {
long nod = 1L;
if (a > b) {
long tmp = a;
a = b;
b = tmp;
}
while (a > 1L && b > 1L) {
for (long i = 2; i <= a; i++) {
if (a % i == 0 && b % i == 0) {
nod *= i;
a /= i;
b /= i;
break;
}
if (a % i == 0) {
a /= i;
break;
}
if (b % i == 0) {
b /= i;
break;
}
}
}
return nod;
}
Первый if отлавливает ситуацию, когда оба числа делятся нацело на переменную цикла и, следовательно, переменная цикла является общим простым (!) множителем для обоих чисел и учитывается для вычисления НОД. Остальные два if отлавливают случаи, когда только одно из чисел делится на переменную цикла; эти множители в НОД не входят.
По-хорошему, цикл for должен перебирать только простые числа. Но нахождение простых чисел является самостоятельной вычислительной задачей, которую здесь решать не хотелось бы. Можно, конечно, использовать таблицу простых чисел, например, в пределах первой тысячи, а для больших чисел, в случае необходимости, вычислять простые числа проверкой на простоту... Но я не стал забираться в дебри факторизации натуральных чисел, а просто закрыл глаза на заведомо холостые проходы цикла for, когда переменная цикла не является простым числом, поскольку после нахождения любого из множителей хотя бы для одного из чисел цикл стартует заново.
Третий вариант по производительности не плох... Но пока это самопал. А что придумали умные люди за многовековую историю математики? "O'key, Гуг... то есть, Википедия, что такое НОД?" Вот, кстати, описано нахождение НОД через каноническое разложение на простые множители, которое мы уже реализовали. А вот что-то новенькое...
04 алгоритм Евклида (рекурсивный)
В самом простом случае алгоритм Евклида применяется к паре положительных целых чисел и формирует новую пару, которая состоит из меньшего числа и разницы между большим и меньшим числом. Процесс повторяется, пока числа не станут равными. Найденное число и есть наибольший общий делитель исходной пары.
Реализуем рекурсивную версию:
long gcd04(long a, long b) {
if (a == b) {
return a;
}
if (a > b) {
long tmp = a;
a = b;
b = tmp;
}
return gcd04(a, b - a);
}
Считается, что рекурсивные алгоритмы менее эффективны, чем итерационные, за счёт накладных расходов на вызов функции. Для проверки делаем и итерационный вариант.
05 алгоритм Евклида (итерационный)
long gcd05(long a, long b) {
while (a != b) {
if (a > b) {
long tmp = a;
a = b;
b = tmp;
}
b = b - a;
}
return a;
}
Кстати, в Викиучебнике есть и другие реализации алгоритма Евклида.
06 бинарный алгоритм (рекурсивный)
НОД(0, n) = n; НОД(m, 0) = m; НОД(m, m) = m;
НОД(1, n) = 1; НОД(m, 1) = 1;
Если m, n чётные, то НОД(m, n) = 2*НОД(m/2, n/2);
Если m чётное, n нечётное, то НОД(m, n) = НОД(m/2, n);
Если n чётное, m нечётное, то НОД(m, n) = НОД(m, n/2);
Если m, n нечётные и n > m, то НОД(m, n) = НОД((n-m)/2, m);
Если m, n нечётные и n < m, то НОД(m, n) = НОД((m-n)/2, n);
Рекурсивная версия:
long gcd06(long a, long b) {
if (a == 0L)
return b;
if (b == 0L)
return a;
if (a == b)
return a;
if (a == 1L || b == 1L)
return 1L;
if (a % 2L == 0L && b % 2L == 0L)
return 2L * gcd06(a / 2L, b / 2L);
if (a % 2L == 0L && b % 2L != 0L)
return gcd06(a / 2L, b);
if (a % 2L != 0L && b % 2L == 0L)
return gcd06(a, b / 2L);
if (a < b)
return gcd06((b - a) / 2L, a);
else
return gcd06((a - b) / 2L, b);
}
07 бинарный алгоритм (итерационный)
Итерационная версия:
long gcd07(long a, long b) {
long nod = 1L;
long tmp;
if (a == 0L)
return b;
if (b == 0L)
return a;
if (a == b)
return a;
if (a == 1L || b == 1L)
return 1L;
while (a != 0 && b != 0) {
if (a % 2L == 0L && b % 2L == 0L) {
nod *= 2L;
a /= 2L;
b /= 2L;
continue;
}
if (a % 2L == 0L && b % 2L != 0L) {
a /= 2L;
continue;
}
if (a % 2L != 0L && b % 2L == 0L) {
b /= 2L;
continue;
}
if (a > b) {
tmp = a;
a = b;
b = tmp;
}
tmp = a;
a = (b - a) / 2L;
b = tmp;
}
if (a == 0)
return nod * b;
else
return nod * a;
}
Кстати, в Викиучебнике есть и другие реализации бинарного алгоритма Евклида.
08 бинарный алгоритм (итерационный) с использованием битовых операций
Поскольку отличительной особенностью бинарного алгоритма является возможность использования битовых сдвигов вместо медленных операций деления и умножения, грех такой возможностью не воспользоваться.
Эта версия функции по логике полностью совпадает с предыдущей, но операции деления и умножения на 2 заменены арифметическими сдвигами, а проверка на чётность – на проверку младшего бита числа.
long gcd08(long a, long b) {
long nod = 1L;
long tmp;
if (a == 0L)
return b;
if (b == 0L)
return a;
if (a == b)
return a;
if (a == 1L || b == 1L)
return 1L;
while (a != 0 && b != 0) {
if (((a & 1L) | (b & 1L)) == 0L) {
nod <<= 1L;
a >>= 1L;
b >>= 1L;
continue;
}
if (((a & 1L) == 0L) && (b & 1L)) {
a >>= 1L;
continue;
}
if ((a & 1L) && ((b & 1L) == 0L)) {
b >>= 1L;
continue;
}
if (a > b) {
tmp = a;
a = b;
b = tmp;
}
tmp = a;
a = (b - a) >> 1L;
b = tmp;
}
if (a == 0)
return nod * b;
else
return nod * a;
}
Ещё немного подготовки
Итак, функции, на которые хватило времени, мозгов и фантазии написаны. Можно допиливать «испытательный стенд» и приступать к тестированию.
Как показали результаты тестовых прогонов, время работы функций сильно зависит от входных данных. Поэтому было решено генерировать случайные пары чисел. Однако выяснилось, что совсем случайные пары чисел слишком часто являются взаимно простыми. Поэтому я решил к паре случайных чисел добавлять случайный множитель. Получилось, что даже если изначально в паре числа были взаимно просты, то алгоритм поиска НОД должен будет вычислить именно этот добавленный множитель.
В остальном алгоритм тестирования прост: для каждого набора случайных данных берётся одна из тестируемых функций и прогоняется REPEAT_TIMES раз. Время работы функции записывается в элемент массива, соответствующий этой функции. После того, как все тесты пройдены, результаты из массива делятся на количество наборов данных – получаем среднее время для набора.
#include <cstdio>
#include <clocale>
#include <cstdlib>
#include <ctime>
using namespace std;
//
// здесь должны быть определения функций вычисления НОД
//
struct sub {
long(*func)(long, long);
const char * comm;
};
sub arr[] = {
{ gcd01, "01 перебор от произвольного числа " },
{ gcd02, "02 перебор от минимального числа " },
{ gcd03, "03 с разложением на делители " },
{ gcd04, "04 алгоритм Евклида рекурсивный " },
{ gcd05, "05 алгоритм Евклида итерационный " },
{ gcd06, "06 бинарный алгоритм рекурсивный " },
{ gcd07, "07 бинарный алгоритм итерационный " },
{ gcd08, "08 бинарный алгоритм итерац. со сдвигом " }
};
const unsigned int RAND_TIMES = 500u;
const unsigned long REPEAT_TIMES = 10000uL;
int main() {
clock_t the_time;
double elapsed_time;
setlocale(LC_ALL, "Russian");
long a, b, c, nod;
srand((unsigned int)time(NULL));
double times[sizeof(arr) / sizeof(sub)] = { 0.0 };
for (unsigned int t = 0u; t < RAND_TIMES; t++) {
c = rand() % 50 + 1;
a = (rand() % 1000 + 1) * c;
b = (rand() % 1000 + 1) * c;
for (unsigned int alg = 0u; alg < sizeof(arr) / sizeof(sub); alg++) {
the_time = clock();
for (unsigned long m = 0uL; m < REPEAT_TIMES; m++) {
nod = arr[alg].func(a, b);
}
elapsed_time = double(clock() - the_time) / CLOCKS_PER_SEC;
times[alg] += elapsed_time;
}
printf("%4u НОД(%7ld, %7ld) = %7ld\n", t + 1, a, b, nod);
}
printf("\nСреднее время для %d пар случайных чисел:\n", RAND_TIMES);
for (unsigned int alg = 0; alg < sizeof(arr) / sizeof(sub); alg++) {
printf("%s: %8.4f сек.\n", arr[alg].comm, times[alg] / RAND_TIMES);
}
return 0;
}
Тестирование
Программа компилировалась под MS Visual Studio 2013 и TDM-GCC 4.8.1 в режиме 64-bit Release.
Среднее время для 500 пар случайных чисел:
-------------------------------------------------------
01 перебор от произвольного числа : 0.5022 сек.
02 перебор от минимального числа : 0.3256 сек.
03 с разложением на делители : 0.0063 сек.
04 алгоритм Евклида рекурсивный : 0.0007 сек.
05 алгоритм Евклида итерационный : 0.0008 сек.
06 бинарный алгоритм рекурсивный : 0.0006 сек.
07 бинарный алгоритм итерационный : 0.0003 сек.
08 бинарный алгоритм итерац. со сдвигом : 0.0002 сек.
Как и ожидалось, первый алгоритм катастрофически неэффективен. Алгоритм №2 – минимальный костыль для №1 – работает почти в 2 раз быстрее.
Третий алгоритм неожиданно показал очень достойный результат: в 50 раз быстрее алгоритма №2.
Четвёртый и пятый варианты – алгоритм Евклида: рекурсивная версия, как ни странно, обогнала итерационную. По сравнению с третьим вариантом время улучшилось почти на порядок.
Бинарный алгоритм Евклида показал наилучшие результаты. Из трёх вариантов реализации рекурсивная версия – самая неторопливая. Наилучший результат у оптимизированной версии с использованием битовых операций.
Итого, самая производительная версия работает более чем в 2500 раз быстрее, чем изначальный вариант.
Примечание. Время, указанное в таблице – это время выполнения REPEAT_TIMES вызовов функции.
Выводы
Выводы достаточно банальны, но всё же:
- Первый пришедший на ум алгоритм решения какой-то задачи часто неоптимален. Хотя и может правильно и достаточно приемлемо работать в определённых условиях.
- Всегда надо попробовать подумать над усовершенствованием алгоритма или над альтернативным вариантом.
- Учите математику!
- Если есть возможность, посмотрите что сделали люди до вас. Может не стоит изобретать велосипед?
- Оптимизируйте код.
Комментарии