интерактивный курс: введение в программирование на PARI/GP

[Dmitriy40]

Aether
Первый цикл - проверяемое число, с шагом 2 чисто для ускорения (чётные числа не подходят по определению).
Присваивание t=x просто инициализирует временную переменную (текущий остаток) для внутренних циклов проверки числа, чтобы не изменить в процессе переменную цикла.
Вектор обнулил для сброса списка делителей.
Второй цикл - проверка делимости на соответствующее число Мерсенна начиная с наибольшего, оно вычисляется и сохраняется в m.
Цикл while - делит (нацело) текущий остаток (в t) на текущее число Мерсенна и расширяет список (вектор) делителей пока текущий остаток ему кратен (t%m==0 вычисляет остаток от деления и проверяет его на 0, а "\" - целочисленное деление).
В итоге t будет равно 1 только если проверяемое число состоит из произведений чисел Мерсенна (в любой положительной степени каждое).
Правда доказать это утверждение я не могу, но как и привести контрпример.

[Aether]

Спасибо, разобрался. Думаю, что код можно ускорить, вычисляя наибольшее число Мерсенна для некоторой области, а не задавая его вручную сразу для всех 100000, так происходит много ненужных сравнений и делений.

[Dmitriy40]

Дык я примерно так и сделал, выставив руками наибольшую степень двойки. Ну а что она больше реально необходимой - не сильно портит, ведь для каждой лишней выполняется фактически одно деление (t%m) и всё, остальное намного быстрее.
Добавление лишнего if(t>=m, ... перед while для пропуска делений если текущий остаток меньше числа Мерсенна ускоряет некритично.
Попробовал предвычислить все числа Мерсенна перед первым циклом - ускорило несущественно.
В общем нет смысла усложнять, значительного ускорения не получается.

[Aether]

Кстати, а как посмотреть и сравнить скорость исполнения кодов?

[Dmitriy40]

Aether в сообщении #1194916 писал(а):

Кстати, а как посмотреть и сравнить скорость исполнения кодов?

После выполнения программы дать команду ## - покажет время выполнения последней команды (или всей последней запущенной программы) в мс.
Правда я выше сравнивал "на глаз".

[Aether]

Что-то 0 мs выдает даже при $10^6$ значений. Как-то неправильно.

[Dmitriy40]

Да, неправильно. Но тогда не знаю. Раньше работал и другой способ, а потом вдруг перестал без причины.
Но вообще с PARI/GP о скорости простых вычислений говорить бессмысленно, у меня счёт до $10^5$ занимал 1.7с, до $10^6$ считалось 21с, а программка на дельфи без оптимизаций за полторы минуты насчитала всё до 2млрд.

[Aether]

У меня новые вопросы. Затеял построить функцию Мёбиуса на натуральном ряду. Для этого необходимо числам, имеющим в факторизации степени простых больше 1 присвоить 0, не имеющим в факторизации степеней простых больше 1, присвоить 1, если количество чисел в факторизации четное и -1, если нечетное. (Прошу не выкладывать готовое решение или выкладывать под спойлер)

Функция factor(n) выдает в качестве результата матрицу из k строк и 2-х столбцов. Чтобы проверить, содержит ли факторизация степени факторов >1, необходимо иметь доступ ко 2-му столбцу матрицы, содержащему степени факторов.

Вопрос, как получить доступ кзначениям 2-го столбца, если матрица задавалась неявно, а формировалась как результат функции factor? Т.е. переменные по которым пробегают индексы матрицы неопределены. Как например извлечь значение из третьей строки второго столбца матрицы?

-- 24.02.2017, 17:50 --

Кажется разобрался, присваиваем переменной m матрицу факторизации: m=factor(n), далее с помощью функции v=matsize(m) определяем её размер, который выдается в виде вектора и который мы сохраняем в переменной v. Первое значение v - количество строк. Далее извлекаем k= v[1] и строим цикл длины k в котором обращаемся к каждой ячейке второго столбца матрицы и сравниваем его с 1: m[k,2]==1. Если все значения 1 - значит факторизация числа не содержит высших степеней, если есть значение больше 1 или !=1, значит содержит высшие степени.

Теперь осталось "самую малость" - воплотить это в работающий код.

[Dmitriy40]

Обрисую интересный вариант проверки вектора на "все единицы" (если меньше 1 значений быть не может) без явных циклов:

v=[...]; if(vecmax(v) == 1, ...)

Получить строку или столбец матрицы в вектор можно примерно так:

v=m[1,]; v=m[,1]

И тогда проверка что число раскладывается на простые множители в степенях не выше 1:

if(vecmax(factor(x)[,2]) == 1, ...)

Красота!

[maxal]

Aether в сообщении #1195040 писал(а):

Затеял построить функцию Мёбиуса на натуральном ряду.

Всё уже построено: moebius(n)

[Aether]

Что-то не заладилось у меня с кодом (

for(n=1,100, m=factor(n); k=matsize(m)[1]; z=0; for(i=1,k, if(m[k,2]==1, z=z+1)); if(z==k, print(n,"=+-1"), print(n,"=0")))

Пропускает некоторые нули

-- 24.02.2017, 19:03 --

maxal в сообщении #1195094 писал(а):

Всё уже построено: moebius(n)

Спасибо, но раз уж начал, то надо доделать, тем более - это важный опыт работы с матрицами и векторами.

-- 24.02.2017, 19:04 --

Dmitriy40 в сообщении #1195093 писал(а):

И тогда проверка что число раскладывается на простые множители в степенях не выше 1:

Ага, спасибо, учту.

[Dmitriy40]

Aether в сообщении #1195101 писал(а):

Что-то не заладилось у меня с кодом (

... if(z==k, print(n,"=+-1)); print(n,"=0"))

Наверное забыли двойную кавычку в первом print.
И второй print выполняется всегда - это так и задумано? Может быть он должен выполняться только если не выполнился первый и тогда следует поместить его в ветку else оператора if?

[Aether]

Все-равно, пропускает некоторые нули, хотя и не выводит значения дважды.

[Dmitriy40]

В if(m[k,2]==1 индекс должен быть не k, а i.

[Aether]

Вот блин, а я уже всю голову сломал, спасибо.