В лекции рассматриваются некоторые арифметические алгоритмы, т.е. алгоритмы над целыми числами.

Будут рассмотрены следующие задачи:

• нахождение НОД;
• возведение числа в натуральную степень;
• проверка числа на простоту;
• факторизация целого числа.

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)

где p₁, p₂, …, p_s – простые числа. Более того, для заданного 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 = 3²×7³×17×19,   

2721411 = 3³×7²×11²×17,

откуда

НОД(997101, 2721411) = 3²×7²×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 – натуральное. Требуется вычислить aⁿ. Прямой способ сделать это – непостредственно перемножить n значений величины a – требует n – 1 операций умножения. Можно ли предложить более эффективный метод? В настоящем разделе мы рассмотрим так называемый бинарный алгоритм возведения в степень, который требует не более 2log2n-1  операций умножения.

Вначале рассмотрим случай, когда n – степень двойки, т. е. n = 2^k для некоторого целого k. Тогда для вычисления aⁿ достаточно выполнить всего k умножений вместо 2^k – 1, требуемых при непосредственном умножении значений a. Пусть, например, n = 16. Формула

a¹⁶ = (((a²)²)²)²

подсказывает нам следующую схему вычислений:

t₁ ← a ∙ a         (t₁ = a²);

t₂ ← t₁ ∙ t₁        (t₂ = a⁴);

t₃ ← t₂ ∙ t₂        (t₃ = a⁸);

t₄ ← t₃ ∙ t₃        (t₄ = a¹⁶).

В случае произвольного n нам поможет его двоичное представление. Пусть, например, n = 19. Поскольку 19 = 10011₂, то a¹⁹ = a^{16 + 2 + 1} = a¹⁶a ²a. Таким образом, надо перемножить значения a¹⁶, a ², a, которые мы уже умеем вычислять. Получаем следующую схему вычисления a¹⁹:

t₁ ← a ∙ a         (t₁ = a²);

t₂ ← t₁ ∙ t₁        (t₂ = a⁴);

t₃ ← t₂ ∙ t₂        (t₃ = a⁸);

t₄ ← t₃ ∙ t₃        (t₄ = a¹⁶);

t₅ ← a ∙ t₁ ∙ t₄    (t₅ = a¹⁹).

Всего мы выполнили 6 операций умножения.

В следующем алгоритме формирование окончательного результата осуществляется параллельно с процедурой получения двоичного представления числа n.

function binary_power_0(a: real, n: longint): real;

{Бинарный алгоритм возведения в степень

(предварительная версия)

Возвращает aⁿ, где n – натуральное}

var p, t: real;

begin

  p := 1.0; {в p формируется результат}

  t := a;   {t принимает значения a, a², a⁴, a⁸, a¹⁶, …}

  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;

{Бинарный алгоритм возведения в степень

Возвращает aⁿ, где n – натуральное}

var p, t: real;

begin

  p := 1.0; {в p формируется результат}

  t := a;   {t принимает значения a, a², a⁴, a⁸, a¹⁶, …}

  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 умножений, однако разложение

a²⁷ = (((a³)³)³

подсказывает, что достаточно 6 умножений (каких?).

Попробуйте предложить схему вычисления a¹⁵, использующую 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, достаточно вычеркивать числа, начиная с p² (а не 2p), так как все составные числа, меньшие p², будут к тому моменту уже вычеркнуты.

В рассматриваемой ниже реализации данного метода мы используем массив логических (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 ).

Здесь p_j называется простым множителем, а k_j – кратностью простого множителя p_j в разложении числа n (обратим внимание, что p_j называется также простым делителем числа 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 с учетом их

               кратностей (всегда достаточно

               длины ≤ log₂ 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.       Темы заданий для самостоятельной работы

• Для всех приведенных функций написать код запускающих программ, предусматривающий ввод и вывод данных. Скомпилировать и запустить программы на ряде примеров.
• Написать рекурсивный вариант бинарного метода вычисления НОД.
• Разработать методы возведения в 15 и 27 степень, требующие соответственно 5 и 6 умножений.
• Разработать и отладить функцию, реализующую метод решета Эратосфена и использующую для представления таблицы чисел тип данных set.

5.       Литература

Арифметическим алгоритмам посвящена большая часть 2-го тома известной монографии Д. Кнута. Помимо рассмотренных в настоящей работе задач и алгоритмов рассматриваются алгоритмы генерации случайных чисел, алгоритмы арифметики произвольной точности, более совершенные алгоритмы проверки простоты и факторизации и др.

Кнут Д. Искусство программирования, том 2. Получисленные методы. – 3-е изд. – М.: Вильямс, 2007. – 832 с.

См. также работу

Ноден П., Китте К. Алгебраическая алгоритмика (с упражнениями и решениями): Пер. с франц. – М.: Мир, 1999. – 720 с.

Некоторые алгоритмы с целыми числами рассматриваются в учебнике

Шень А. Программирование: теоремы и задачи. – М.: МЦНМО, 1995. – 264 с.

Вопросам низкоуровневой оптимизации алгоритмов с целыми числами посвящена работа

Уоррен Г. С. Алгоритмические трюки для программистов: Пер. с англ. – М.: Вильямс, 2003. – 288 с.