Отправляет email-рассылки с помощью сервиса Sendsay
  Все выпуски  

Система компьютерной алгебры GAP - # 29, 28-10-03


Информационный Канал Subscribe.Ru

Рассылка "Система компьютерной алгебры GAP"
ведущий рассылки А.Б.Коновалов, a_konovalov@hotmail.com
выпуск 29 от 28 октября 2003 г.

На сайте Украинской группы пользователей GAP в разделе "Изучаем алгебру с GAP" опубликован материал посвященный теории чисел: элементарным функциям, имеющимся в системе компьютерной алгебры GAP для работы с целыми числами, алгоритму Евклида, простым числам и разложению целых чисел на множители. При его подготовке использована дипломная работа студентки Запорожского государственного университета Н.В.Жатько (2003 г.). Полный текст приводится ниже.

Некоторые функции для работы с целыми числами

Функция IsInt(obj) возвращает true или false в зависимости от того, является ли obj целым числом, или нет. Аналогично, IsPosInt(obj) проверяет, является ли obj целым положительным числом. Функция Int(elm) возвращает целое число int, чье значение зависит от типа входного параметра elm. Если elm – рациональное число, то int – его целая часть. Если elm – элемент конечного поля простой характеристики, то int – наименьшее неотрицательное число, такое что elm = int * One(elm). Если elm – строка, состоящая из цифр 0, 1, ..., 9 и, возможно, знака "минус" в первой позиции, то int – описываемое этой строкой целое число, например:
gap> Int(4/3); Int(-2/3);
1
0
gap> int:= Int(Z(5)); int * One(Z(5));
2
Z(5)
gap> Int("12345"); Int("-27"); Int("-27/3");
12345
-27
fail

Функция IsEvenInt(n) проверяет, является ли число четным. Аналогично, функция IsOddInt(n) проверяет, является ли число нечетным. Функция AbsInt(n) возвращает модуль целого числа n:
gap> AbsInt(33);
33
gap> AbsInt(-214378);
214378
gap> AbsInt(0);
0

Функция SignInt(n) возвращает знак числа n.
gap> SignInt(33);
1
gap> SignInt(-214378);
-1
gap> SignInt(0);
0

Функция LogInt(n, base) возвращает целую часть логарифма положительного целого числа n по положительному целому основанию base. Если n или base не положительны, будет получено сообщение об ошибке.
gap> LogInt(1030, 2);
10 # так как 2^10 = 1024
gap> LogInt(1, 10);
0

Функции RootInt(n) и RootInt(n, k) возвращают целую часть корня k-й степени из целого числа n. Если k не задано, то вычисляется квадратный корень. Если n положительно, то RootInt возвращает наибольшее положительное целое r, такое что r/k меньше или равняется n. Если n отрицательно и k нечетно, то RootInt возвращает -RootInt(-n,k). Если n отрицательно и k четно, то возвращается сообщение об ошибке. Оно также будет выдано, если k равно нулю или отрицательно.
gap> RootInt(361);
19
gap> RootInt(2 * 10^12);
1414213
gap> RootInt(17000, 5);
7 # так как 7^5 = 16807

Функция SmallestRootInt(n) возвращает наименьший корень из целого числа n, т.е. минимальное по модулю целое число r, для которого существует положительное число k, такое что n = r^k.
gap> SmallestRootInt(2^30);
2
gap> SmallestRootInt(-(2^30));
-4 # так как (-2)^30 = 2^30
gap> SmallestRootInt(279936);
6
gap> LogInt(279936, 6);
7
gap> SmallestRootInt(1001);
1001

Функция Random(Integers) возвращает псевдослучайное целое число в диапазоне от -10 до 10. Если необходимы случайные числа из другого интервала, используется функция Random([low .. high]).


Алгоритм Евклида

Вначале заметим, что в GAP имеются стандартные функции Gcd и Gcdex для нахождения наибольшего общего делителя, а также линейного представления НОД и нуля:
gap> Gcd(100, 48);
4
gap> Gcdex(100, 48);
rec( gcd := 4, coeff1 := 1, coeff2 := -2, coeff3 := -12, coeff4 := 25 )
gap> 100*1+48*-2;
4
gap> 100*-12+48*25;
0
Однако, использование этих функций выводит на экран только конечный результат. Поэтому покажем, как, например, разработать программу, с помощью которой мы сможем получить информацию о каждом этапе алгоритма Евклида.

Вначале разработаем функцию для вычисления неполного частного от деления a на b, которая будет находить остаток от деления a на b, вычитать его из a, а затем делить полученное число на b (заметим, что стандартная функция системы QuoInt возвращает неполное частное при делении модуля a на модуль b со знаком, равным произведению знаков a и b, так что QuoInt(-5,3) = -1, тогда как quotient(-5,3) = -2):
quotient:=function(a,b)
local r, m, q;
# a = b * q + r
r:= a mod b; # остаток от деления
m:= a - r; # m = b * q
q:=m/b;
return q;
end;

Пример работы с данной функцией:
gap> quotient(18,3);
6
gap> quotient(27,8);
3
gap> quotient(2^100-1,2^100+1);
0
gap> quotient(2^100+1,2^100-1);
1
Теперь реализуем алгоритм Евклида. Строго говоря, предыдущая функция не является необходимой для него, так как на каждом шаге нас будет интересовать только остаток:
euclid:=function(a,b)
local r;
r:=a mod b;
while r <> 0 do
a:=b;
b:=r;
r:=a mod b;
od;
return b;
end;

Пример вычисления НОД с помощью данной функции:
gap> euclid(21,9);
3
gap> euclid(7,9);
1
gap> euclid(2^100-1,2^100+1);
1
gap> euclid(6456,46);
2
Однако, для того, чтобы поэтапно вывести на печать алгоритм Евклида, нам все же понадобится дополнительное обращение к разработанной ранее функции quotient. Для этого модифицируем функцию euclid следующим образом:
euclidPrint:=function(a,b)
local r,q;
r:=a mod b;
q:=quotient(a,b);
Print(a, " = ", b, " * ", q, " + ", r, "\n");
while r <> 0 do
a:=b;
b:=r;
r:=a mod b;
q:=quotient(a,b);
Print(a, " = ", b, " * ", q, " + ", r, "\n");
od;
Print("NOD = ", b, "\n");
return b;
end;

Теперь алгоритм Евклида будет выводиться на печать:
gap> euclidPrint(300,85);
300 = 85 * 3 + 45
85 = 45 * 1 + 40
45 = 40 * 1 + 5
40 = 5 * 8 + 0
NOD = 5
5
gap> euclidPrint(4564300,6545685);
4564300 = 6545685 * 0 + 4564300
6545685 = 4564300 * 1 + 1981385
4564300 = 1981385 * 2 + 601530
1981385 = 601530 * 3 + 176795
601530 = 176795 * 3 + 71145
176795 = 71145 * 2 + 34505
71145 = 34505 * 2 + 2135
34505 = 2135 * 16 + 345
2135 = 345 * 6 + 65
345 = 65 * 5 + 20
65 = 20 * 3 + 5
20 = 5 * 4 + 0
NOD = 5
5


Простые числа и разложение на множители

С помощью функции Primes может быть получен список всех простых чисел, не превышающих 1000:
gap> Primes;
[2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997]

Он содержит 168 чисел:
gap> Length(last);
168

Так как Primes – это список, то мы можем, например, получить 20-е простое число следующим образом:
gap> Primes[20];
71

Для того, чтобы проверить, является ли число простым, используется функция IsPrimeInt (этот же результат можно получить и с помощью более общей функции IsPrime, которая применима не только к целым числам, но и к элементам Евклидовых колец):
gap> IsPrimeInt(657867897897897897891);
false
gap> 2^31-1;
2147483647
gap> IsPrimeInt(2^31-1);
true
gap> IsPrimeInt(2^256-1);
false

0 и 1 не являются простыми числами, и GAP "знает" это:
gap> IsPrimeInt(0);
false
gap> IsPrimeInt(1);
false
Здесь и далее функции могут работать и с отрицательными числами:
gap> IsPrimeInt(-7);
true
gap> IsPrimeInt(-10);
false

Для чисел, больших чем 1013, используемый алгоритм проверки простоты числа может выдавать подобное предупреждение:
gap> a:=657867897897321867697897897891;
657867897897321867697897897891
gap> IsPrimeInt(a);
#I beyond the guaranteed bound of the probabilistic primality test
true

Однако до сих пор неизвестны составные числа, большие чем 1013, которые этим тестом определяются как простые. Вывод этого предупреждения можно отключить, если использовать функцию IsProbablyPrimeInt:
gap> IsProbablyPrimeInt(a);
true

Функция IsPrimePowerInt определяет, является ли число степенью простого числа:
gap> IsPrimePowerInt(1000);
false
gap> IsPrimePowerInt(2^31);
true
gap> IsPrimePowerInt(2^31-1);
true

Заметим, что 231-1 – простое число (см. выше), что и объясняет последний результат.

Функции NextPrimeInt и PrevPrimeInt возвращают соответственно минимальное простое число, которое больше заданного, и максимальное простое число, которое меньше заданного.
gap> NextPrimeInt(2^32);
4294967311
gap> PrevPrimeInt(2^32);
4294967291
gap> NextPrimeInt(2^32)-2^32;
15
gap> 2^32-PrevPrimeInt(2^32);
5

Для разложения на множители используется функция FactorsInt (этот же результат можно получить и с помощью более общей функции Factors, которая применима не только к целым числам, но и к элементам Евклидовых колец):
gap> FactorsInt(2^64-1);
[3, 5, 17, 257, 641, 65537, 6700417]
gap> FactorsInt(2^128-1);
[3, 5, 17, 257, 641, 65537, 274177, 6700417, 67280421310721]
gap> FactorsInt(2^200-1);
[3, 5, 5, 5, 11, 17, 31, 41, 101, 251, 401, 601, 1801, 4051, 8101, 61681, 268501, 340801, 2787601, 3173389601]

К этому разложению можно применить функцию Collected для его сокращения:
gap> a:=2^20*3^15*5^7;
1175462461440000000
gap> FactorsInt(a);
[2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5]
gap> Collected(FactorsInt(a));
[[2, 20], [3, 15], [5, 7]]

Можно использовать и функцию PrimePowersInt для представления этого списка в другом виде:
gap> PrimePowersInt(a);
[2, 20, 3, 15, 5, 7]

Для вывода на печать канонического разложения составного числа на множители используется функция PrintFactorsInt:
gap> PrintFactorsInt(a);
2^20*3^15*5^7
gap> PrintFactorsInt(Factorial(20));
2^18*3^8*5^4*7^2*11*13*17*19

Функция DivisorsInt возвращает список всех положительных делителей целого числа:
gap> DivisorsInt(100);
[1, 2, 4, 5, 10, 20, 25, 50, 100]

Функции Sigma и Tau вычисляют соответственно сумму и число всех положительных делителей целого числа:
gap> Sigma(100);
217
gap> Tau(100);
9
gap> Length(DivisorsInt(100))=Tau(100);
true
gap> Sum(DivisorsInt(100))=Sigma(100);
true

С уважением,
Коновалов Александр Борисович , председатель Украинской группы пользователей GAP ,
доцент кафедры алгебры и геометрии Запорожского государственного университета




http://subscribe.ru/
E-mail: ask@subscribe.ru
Отписаться
Убрать рекламу

В избранное