|
В лекции рассматриваются некоторые арифметические алгоритмы, т.е. алгоритмы над целыми числами. Будут рассмотрены следующие задачи:
1. Наибольший общий делительВ настоящем разделе мы рассмотрим задачу нахождения наибольшего общего делителя двух целых чисел и познакомимся с некоторыми алгоритмами решения этой задачи. 1.1. ОпределенияГоворят, что целое число d является делителем целого числа u, если для некоторого целого q справедливо u = qd. В этом случае также говорят, что u кратно d. Например, 3, −5, 30 являются делителями числа 30, а 7 таковым не является. Если некоторое целое число d является делителем одновременно для каждого из двух заданных чисел u и v, то d называется их общим делителем. Например, общими делителями чисел 12 и 30 являются 1, 2, 6, −3 и др., а, скажем, 5 их общим делителем не является. Среди всех общих делителей двух заданных целых чисел u и v выберем максимальный. Он называется наибольшим общим делителем (сокращенно НОД) этих чисел и обозначается НОД(u, v) или gcd(u, v) (сокращение gcd – от greatest common divisor). Например, НОД(12, 30) = 6. Очевидно, НОД определен, если по крайней мере одно из чисел u и v не равно нулю. Обычно полагают НОД(0, 0) = 0. Легко видеть, что НОД(u, v) = НОД(|u|, |v|) и НОД(u, 0) = |u|, поэтому достаточно рассматривать задачу нахождения НОД для натуральных u и v. 1.2. Нахождение НОД с помощью разложения на множителиТрадиционно изучаемый в школе метод нахождения НОД основан на использовании так называемого канонического разложения целого числа. Напомним, что целое число p, большее 1, называется простым, если его делителями являются только ±1, ±p. Следующие числа являются простыми: 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, … Простых чисел бесконечно много. Целое число, большее 1, не являющееся простым, называется составным. Целые числа 0, 1 и отрицательные целын числа мы не будет относить ни к простым, ни к составным. Хорошо известно (основная теорема арифметики), что любое натуральное число u допускает представление в виде u=p1k1×p2k2×…×psks , p1<p2<…<ps (1) где p1, p2, …, ps – простые числа. Более того, для заданного u представление (1) единственно. Представление (1) называется разложением числа u на простые множители или каноническим разложением. Если канонические разложения чисел u и v известны, то для нахождения НОД(u, v) достаточно отобрать общие простые множители: если канонические разложения u и v суть u=p1k1×p2k2×…×psks , v=p1l1×p2l2×…×psls , тогда НОД(u, v)=p1min{k1,l1}×p2min{k2,l2}×…×psmin{ks,ls} . Пусть, например, требуется найти НОД(997101, 2721411). Находим их канонические разложения: 997101 = 32×73×17×19, 2721411 = 33×72×112×17, откуда НОД(997101, 2721411) = 32×72×17 = 7497. Очевидный недостаток данного метода заключается в его крайней медлительности. Наиболее трудоемкий процесс здесь – поиск канонического разложения. На самом деле, задача нахождения канонического разложения целого числа (называемая также задачей факторизации) является намного более трудной, чем задача нахождения НОД. Сведение задачи к более трудоемкой нецелесообразно. Задача факторизации рассматривается в разделе 3.3. 1.3. Алгоритм ЕвклидаВ настоящем разделе мы познакомимся с намного более производительным, чем метод из раздела 1.2, алгоритмом нахождения НОД – алгоритмом Евклида, но прежде напомним еще некоторые свойства целых чисел. Пусть u – целое, а v – натуральное число. Тогда существуют, причем единственные, целые числа q и r, такие, что u = qv + r, 0 ≤ r < v. (2) Правило, которое произвольной паре целых чисел u, v, где v ≠ 0, ставит в соответствие пару целых чисел q и r, удовлетворяющих соотношениям (2), называется операцией деления с остатком, причем u называется делимым, v – делителем, q – неполным частным, а r – остатком. Для неполного частного q и остатка r будем использовать обозначения q = u div v, r = u mod v. Приведем несколько примеров: 14 div 5 = 2, 14 mod 5 = 4, так как 14 = 2×5 + 4, 0 ≤ 4 < 5; (−14) div 5 = −3, (−14) mod 5 = 1, так как −14 = (−3)×5 + 1, 0 ≤ 1 < 5; 3 div 7 = 0, 3 mod 7 = 3, так как 3 = 0×7 + 3, 0 ≤ 3 < 7; (−3) div 7 = −1, (−3) mod 7 = 4, так как −3= (−1)×7 + 4, 0 ≤ 4 < 7. Алгоритм Евклида основан на последовательном использовании легко доказываемого соотношения НОД(u, v) = НОД(v, u mod v). (3) Пусть, например, требуется найти НОД(88179, 5313). Применяя (3) многократно, получаем НОД(88179, 5313) = = НОД(5313, 88179 mod 5313) = НОД(5313, 3171) = = НОД(3171, 5313 mod 3171) = НОД(3171, 2142) = = НОД(2142, 3171 mod 2142) = НОД(2142, 1029) = = НОД(1029, 2142 mod 1029) = НОД(1029, 84) = = НОД(84, 1029 mod 84) = НОД(84, 21) = = НОД(21, 84 mod 21) = НОД(21, 0) = 21. Мы приходим к следующему алгоритму. Вначале запишем рекурсивный вариант, который, по большему счету, представляет непосредственную запись формулы (3) и равенства НОД(u, 0) = u: function euclid_recursive(u, v) {Алгоритм Евклида (рекурсивный вариант) Возвращает НОД(u, v), где u, v – целые, u ≥ 0, v ≥ 0} begin if v = 0 do euclid_recursive = u; else euclid_recursive = euclid_recursive(v, u mod v); end;
Нетрудно записать и итеративный алгоритм: function euclid_0(u, v: longint): longint; {Алгоритм Евклида (предварительный вариант) Возвращает НОД(u, v), где u, v – целые, u ≥ 0, v ≥ 0} var tmp: longint; begin while v > 0 do begin tmp := v; v := u mod v; u := tmp; end; euclid_0 := u; end;
Добавим две строчки, чтобы метод работал для любых целых (а не только неотрицательных): function euclid(u, v: longint): longint; {Алгоритм Евклида Возвращает НОД(u, v), где u, v – целые} var tmp: longint; begin u = abs(u); v = abs(v); while v > 0 do begin tmp := v; v := u mod v; u := tmp; end; euclid := u; end;
1.4. Бинарный алгоритмЕще один известный метод вычисления НОД – так называемый бинарный алгоритм. Он основан на следующих простых фактах: если u и v – четные, то НОД(u, v) = 2 НОД(u/2, v/2); если u – четное, v – нечетное, то НОД(u, v) = НОД(u/2, v); если u и v – нечетные, то НОД(u, v) = НОД(u − v, v) = НОД((u – v)/2, v). function binary_gcd_0(u, v: longint): longint; {Бинарный алгоритм нахождения НОД (предварительная версия) Возвращает НОД(u, v), где u, v – целые, u ≥ 0, v ≥ 0} var d: longint; begin d := 1; while (u > 0) and (v > 0) do begin if u mod 2 = 0 then if v mod 2 = 0 then begin u := u div 2; v := v div 2; d := 2*d; end else u := u div 2; else if v mod 2 = 0 then v := v div 2; else if u > v then u := (u – v) div 2; else v := (v – u) div 2; end; if u = 0 then binary_gcd_0 = d*v; else binary_gcd_0 = d*u; end; Заметим, что если на некоторой итерации алгоритма одно из чисел u и v нечетно, то на следующих итерациях u и v уже не могут стать одновременно четными. Это позволяет вынести вычисление множителя d за пределы основного цикла. Мы приходим к следующему усовершенствованному алгоритму. function binary_gcd(u, v: longint): longint; {Бинарный алгоритм нахождения НОД Возвращает НОД(u, v), где u, v – целые, u ≥ 0, v ≥ 0} var d: longint; begin if (u = 0) and (v = 0) then binary_gcd = 0; else begin d := 1; while (u mod 2 = 0) and (v mod 2 = 0) do begin u := u div 2; v := v div 2; d := 2*d; end while (u > 0) and (v > 0) do begin if u mod 2 = 0 then u := u div 2; else if v mod 2 = 0 then v := v div 2; else if u > v then u := (u – v) div 2; else v := (v – u) div 2; end; if u = 0 then binary_gcd = d*v; else binary_gcd = d*u; end; end;
Мы предоставляем читателю записать рекурсивный вариант бинарного алгоритма нахождения НОД. 2. Бинарный алгоритм возведения в степеньРассмотрим задачу возведения действительного числа в натуральную степень. Пусть a – действительное число, а n – натуральное. Требуется вычислить an. Прямой способ сделать это – непостредственно перемножить n значений величины a – требует n – 1 операций умножения. Можно ли предложить более эффективный метод? В настоящем разделе мы рассмотрим так называемый бинарный алгоритм возведения в степень, который требует не более 2log2n-1 операций умножения. Вначале рассмотрим случай, когда n – степень двойки, т. е. n = 2k для некоторого целого k. Тогда для вычисления an достаточно выполнить всего k умножений вместо 2k – 1, требуемых при непосредственном умножении значений a. Пусть, например, n = 16. Формула a16 = (((a2)2)2)2 подсказывает нам следующую схему вычислений: t1 ← a ∙ a (t1 = a2); t2 ← t1 ∙ t1 (t2 = a4); t3 ← t2 ∙ t2 (t3 = a8); t4 ← t3 ∙ t3 (t4 = a16). В случае произвольного n нам поможет его двоичное представление. Пусть, например, n = 19. Поскольку 19 = 100112, то a19 = a16 + 2 + 1 = a16a 2a. Таким образом, надо перемножить значения a16, a 2, a, которые мы уже умеем вычислять. Получаем следующую схему вычисления a19: t1 ← a ∙ a (t1 = a2); t2 ← t1 ∙ t1 (t2 = a4); t3 ← t2 ∙ t2 (t3 = a8); t4 ← t3 ∙ t3 (t4 = a16); t5 ← a ∙ t1 ∙ t4 (t5 = a19). Всего мы выполнили 6 операций умножения. В следующем алгоритме формирование окончательного результата осуществляется параллельно с процедурой получения двоичного представления числа n. function binary_power_0(a: real, n: longint): real; {Бинарный алгоритм возведения в степень (предварительная версия) Возвращает an, где n – натуральное} var p, t: real; begin p := 1.0; {в p формируется результат} t := a; {t принимает значения a, a2, a4, a8, a16, …} while n > 0 do begin if n mod 2 = 1 then p := p*t; {если n нечетно, то умножаем p на t} n := n div 2; {переходим к следующему биту} t := t*t; end; binary_power_0 = p; end; На последней итерации алгоритм выполняет одно «лишнее» умножение t := t*t (найденное значение t далее не используется). Устраним этот недостаток. function binary_power(a: real, n: longint): real; {Бинарный алгоритм возведения в степень Возвращает an, где n – натуральное} var p, t: real; begin p := 1.0; {в p формируется результат} t := a; {t принимает значения a, a2, a4, a8, a16, …} while true do begin if n mod 2 = 1 then p := p*t; {если n нечетно, то умножаем p на t} n := n div 2; {переходим к следующему биту} if n = 0 then break; t := t*t; end; binary_power = p; end; Количество итераций данного алгоритм равно числу разрядов в двоичной записи числа n, т. е. log2n . На каждой итерации выполняется не более 2 умножений, а напоследней итерации – одно. Получаем, что общее число умножений в бинарном алгоритме возведения в степень не превосходит 2log2n-1 ; точное значение равно log2n+nbits(n)-1 , где nbits(n) – количество ненулевых разрядов в двоичном представлении числа n. Это значительно меньше количества умножений, n – 1, требуемых при непосредственном умножении значений a. Несмотря на то, что метод весьма быстрый, не для всякого n он использует минимально возможное количество умножений. Так, например, для n = 27 бинарный метод использует 7 умножений, однако разложение a27 = (((a3)3)3 подсказывает, что достаточно 6 умножений (каких?). Попробуйте предложить схему вычисления a15, использующую 5 умножений, а не 6 умножений, как в бинарном алгоритме. 3. Проверка на простоту. Факторизация3.1. Проверка целого числа на простотуПусть задано натуральное число n. Требуется определить, простое ли оно. Тривиальный алгоритм, вытекающий их определения простого числа, состоит в переборе чисел от 2 до n – 1 и поиске среди них делителей числа n. Если делитель найден, то стоп – число n составное. Если делителя среди всех рассмотренных чисел нет, то n простое. Процедура использует не более n – 2 пробных делений. Конечно, этот метод можно улучшить! Во-первых, в качестве пробных делителей достаточно рассматривать числа, не превышающие n . Действительно, если у n есть делитель k>n , то его делителем будет также число n/k≤n . Во-вторых, при проверке можно опускать четные числа, большие двух, так как все они составные (если число n составное, то у него всегда есть простой делитель!). Мы приходим к следующему алгоритму пробных делений: function is_prime_0(n: longint): boolean; { Проверка числа на простоту (предварительный вариант) Вход: n - заданное целое число, n > 1 Функция возвращает true, если n - простое false, если n - составное } var d, dmax: longint; begin if n = 2 is_prime_0 := true; else if n mod 2 = 0 is_prime_0 := false; else begin is_prime_0 := true; dmax := sqrt(n); d := 3; {d – пробный делитель} while d < dmax do begin if n mod d = 0 then begin is_prime_0 := false; break; end d = d + 2; {следующее нечетное число} end; end; Данный алгоритм использует не более 12n пробных делений. Трудоемкость можно сократить еще приблизительно в полтора раза и довести до (приблизительно) 13n пробных делений. Для этого заметим, что все простые числа, кроме 2 и 3, имеют вид 6k ± 1, где k – натуральное, поэтому в качестве пробных делителей достаточно брать только такие числа. Мы приходим к следующему усовершенствованному алгоритму пробных делений: function is_prime(n: longint): boolean; { Проверка числа на простоту Вход: n - заданное целое число, n > 1 Функция возвращает true, если n - простое false, если n - составное } var d, k, kmax: longint; if (n = 2) or (n = 3) or (n = 5) or (n = 7) is_prime := true; else if n mod 2 = 0 is_prime := false; {n четное} else if n mod 3 = 0 is_prime := false; {n кратно 3} else begin is_prime := true; kmax := (sqrt(n) + 1) div 6; for k := 1 to kmax do begin d := 6*k - 1; {пробный делитель} if n mod d = 0 begin is_prime := false; break; end; d := 6*k + 1; {пробный делитель} if n mod d = 0 begin is_prime := false; break; end; end; end; end;
3.2. Таблица простых чисел. Решето ЭратосфенаТеперь рассмотрим задачу построения списка простых чисел. Пусть требуется найти k первых простых числа. Мы знаем, что 2 – простое. Далее рассматриваем все нечетные числа p в порядке возрастания и проверяем каждое из них на простоту. Для проверки будем действовать, как и в алгоритме пробных делений, но теперь в качестве пробных делителей будем рассматривать только найденные к настоящему моменту простые числа, не превышающие p . Итерации повторяем, пока не найдем k простых чисел. procedure primes(k: longint, var primes: array of longint) { Таблица простых чисел.
Функция находит k первых простых чисел. Вход: k - количество искомых простых чисел, k ≥ 2 primes – массив длины не менее k Выход: primes содержит k первых простых чисел } var p, i, j, pj, sqrtp: longint; next_prime_found: boolean; begin primes[1] = 2; primes[2] = 3; p = 3; {последнее найденное простое число} for i := 3 to k do repeat {ищем следующее простое число} p := p + 2; {следующее испытуемое число} sqrtp := sqrt(p); j := 1; while true do {перебираем пробные делители} begin pj := primes[j]; {пробный делитель} if p mod pj = 0 then begin next_prime_found := false; break; end; if pj >= sqrtp then begin primes[j] := p; {p - простое} next_prime_found := true; break; end; j = j + 1; end; until next_prime_found; end;
Познакомимся с другим методом построения списка простых чисел, называемым решетом Эратосфена. Вначале рассмотрим пример. Предположим, что нам требуется найти все простые числа, не превосходящие 30. Выпишем все числа от 2 до 30: 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Мы знаем, что 2 – простое число. Вычеркнем из списка все четные числа (кроме 2): 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Следующее невычеркнутое число – это 3. Оно простое. Вычеркиваем все большие числа, кратные 3 (некоторые из них уже были вычеркнуты ранее): 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 В полученном списке находим следующее невычеркнутое число – это 5. Оно простое. Вычеркиваем все большие числа, кратные 5 (некоторые из них уже были вычеркнуты ранее): 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 Мы утверждаем, что невычеркнутыми останутся все простые числа и только они (почему?). Рассмотрим общий случай: требуется найти все простые числа, не превосходящие n. Выписываем все числа от 1 до n и начинаем выполнять аналогичные действия: вначале вычеркиваем все числа кратные и превосходящие 2, затем кратные и превосходящие 3 и т. д. На каждой итерации необходимо найти следующее невычеркнутое число p: оно обязятельно будет простым, а затем вычеркнуть все кратные ему (и превосходящие его) числа, т. е. 2p, 3p, 4p и т. д. до конца списка. Итерации достаточно продолжать, пока p≤n . По окончании итераций список будет содержать все простые числа, не большие n, и только их. Заметим, что когда рассматривается число p, достаточно вычеркивать числа, начиная с p2 (а не 2p), так как все составные числа, меньшие p2, будут к тому моменту уже вычеркнуты. В рассматриваемой ниже реализации данного метода мы используем массив логических (boolean) значений b[1..n], причем b[j] = true тогда и только тогда, когда число j не вычеркнуто из списка (элемент b[1] не используется). procedure Eratosthenes(n: integer, var b: array of boolean) { Решето Эратосфена. Вход: n – целое число, n ≥ 1 b – логический массив длины не менее n Выход: b содержит логические значения, такие, что b[j] = true тогда и только тогда, когда j - простое } var p, j, sqrtn: int; begin b[1] := false; for j := 2 to n do {инициализация массива} b[j] := true; {ни одно число не вычеркнуто} sqrtn := sqrt(n); for p := 2 to sqrtn do if b[p] then {если p не вычеркнуто, т.е. p - простое} for j := p*p to n step p do b[j] := false; {вычеркиваем числа, кратные p} end;
3.3. Факторизация целых чиселВ настоящем разделе рассмотрим задачу факторизации целого числа, т. е. задачу разложения заданного целого числа n ≥ 2 в виде произведения простых чисел n=p1k1×p2k2×…×psks , p1<p2<…<ps , kj≥1 (j=1, 2,…, s ). Здесь pj называется простым множителем, а kj – кратностью простого множителя pj в разложении числа n (обратим внимание, что pj называется также простым делителем числа n). Рассмотрим простой алгоритм факторизации. Пусть задано число n ≥ 2. С помощью пробных делений ищем его наименьший простой делитель p. Если p < n, то делим n на p и повторяем поиск простых делителей для нового n. Заметим, что поиск делителя для нового n достаточно возобновить c последнего найденного p. Итерации заканчиваются, как только устанавливается, что n простое. В приведенной ниже процедуре в качестве пробных делителей будем рассматривать числа 2, 3 и числа вида 6k ± 1, где k – натуральное. procedure factor(n: longint, var divisors: array of longint, var nd: longint) { Разложение числа на простые множители.
Функция находит все простые делители заданного числа.
Вход: n - заданное число, n ≥ 2 divisors – массив длины, достаточной для хранения всех простых делителей числа n с учетом их кратностей (всегда достаточно длины ≤ log2 n) Выход: divisors содержит простые делители числа n; каждый делитель повторяется столько раз, какова его кратность; делители расположены в порядке возрастания nd - количество делителей с учетом их кратностей } var p, k: longint; begin nd := 0; while n mod 2 = 0 do begin nd := nd + 1; divisors[nd] := 2; n := n div 2; end while n mod 3 = 0 do begin nd := nd + 1; divisors[nd] := 3; n := n div 3; end k := 1; while true do begin p := 6*k - 1; {пробный делитель} if n mod p = 0 then begin nd := nd + 1; divisors[nd] := p; n := n div p; end; if p*p > n then break; p := 6*k + 1; {пробный делитель} if n mod p = 0 then begin nd := nd + 1; divisors[nd] := p; n := n div p; end; if p*p > n then break; end; nd := nd + 1; divisors[nd] := n; end;
4. Темы заданий для самостоятельной работы
5. ЛитератураАрифметическим алгоритмам посвящена большая часть 2-го тома известной монографии Д. Кнута. Помимо рассмотренных в настоящей работе задач и алгоритмов рассматриваются алгоритмы генерации случайных чисел, алгоритмы арифметики произвольной точности, более совершенные алгоритмы проверки простоты и факторизации и др. Кнут Д. Искусство программирования, том 2. Получисленные методы. – 3-е изд. – М.: Вильямс, 2007. – 832 с. См. также работу Ноден П., Китте К. Алгебраическая алгоритмика (с упражнениями и решениями): Пер. с франц. – М.: Мир, 1999. – 720 с. Некоторые алгоритмы с целыми числами рассматриваются в учебнике Шень А. Программирование: теоремы и задачи. – М.: МЦНМО, 1995. – 264 с. Вопросам низкоуровневой оптимизации алгоритмов с целыми числами посвящена работа Уоррен Г. С. Алгоритмические трюки для программистов: Пер. с англ. – М.: Вильямс, 2003. – 288 с.
|
|
© a-urusov2011 |