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

[Dmitriy40]

gris в сообщении #1626336 писал(а):

Никак не мог быстро находить максимальный простой делитель не очень больших чисел.

Прямой функции нет. Это вообще-то по сложности эквивалентно разложению на множители, что как известно сложная задача.
Но улучшить код можно: начинать проверку от корня из числа, лучше округлённого вниз до нечётного числа и идти с шагом -2 (чётность проверить сразу отдельно).
А можно сохранить все простые до сотен тысяч в массив (десяток тысяч чисел) командой primes() и идти по нему с конца (а точку начала обхода определить опять же как корень из числа, для сортированных массивов setsearch работает быстро и с ненулевым третьим аргументом выдаст и точку отсутствующего значения, что и надо).

gris в сообщении #1626336 писал(а):

И ещё: не грешно ли изменять внутри цикла переменную цикла? Например, вместо break или next(3)?

Некоотрые циклы это позволяют (for, forstep), некоторые не позволяют (forprime). Но это очень плохая идея! Лучше никогда так не делать. Структурное программирование придумали не от балды и не от хорошей жизни, это работающая методика уменьшения человеческих ошибок при написании кода.
Хочется менять переменную внутри цикла - используйте while или даже любой бесконечный цикл (например while(1,... ) и делайте внутри с переменными что хотите.

[wrest]

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

Но улучшить код можно: начинать проверку от корня из числа, лучше округлённого вниз до нечётного числа и идти с шагом -2

Хм... максимальный простой делитель может быть больше чем корень...

[Dmitriy40]

wrest в сообщении #1626347 писал(а):

Ну я бы тогда начинал с простого, не большего корня n.

Несущественно.

wrest в сообщении #1626347 писал(а):

и шёл бы вниз по опять же precprime() а не просто по нечётным.

Вы правы, пока простых хватает, то так быстрее:

? default(primelimit)
%1 = 51000000
? x=5*3^30;forstep(p=precprime(sqrtint(x)),1,-2,if(x%p==0&&ispseudoprime(p),print(x/1e15,"e15:",p);break));
1.0294556604732450000000000000000000000e15:5
time = 1,826 ms.
? x=5*3^30;p=sqrtint(x)+1;while((p=precprime(p-1))>2,if(x%p==0,print(x/1e15,"e15:",p);break));
1.0294556604732450000000000000000000000e15:5
time = 3,230 ms.

Но когда простых перестаёт хватать, то ситуация быстро меняется на обратную:

? x=53*3^30;forstep(p=precprime(sqrtint(x)),1,-2,if(x%p==0&&ispseudoprime(p),print(x/1e15,"e15:",p);break));
10.912230001016397000000000000000000000e15:53
time = 6,006 ms.
? x=53*3^30;p=sqrtint(x)+1;while((p=precprime(p-1))>2,if(x%p==0,print(x/1e15,"e15:",p);break));
10.912230001016397000000000000000000000e15:53
time = 10,858 ms.

wrest в сообщении #1626348 писал(а):

Хм... максимальный простой делитель может быть больше чем корень...

Хм, это я походу с минимальным спутал ... Ну тогда не с корня, а с половины числа.
Но соотношения почти не изменятся. только порог чисел станет не квадрат лимита, а сам удвоенный лимит.

-- 18.01.2024, 02:39 --

gris
Вот только для таких малых чисел, у которых именно минимальный простой делитель не более сотен тысяч, конструкция f=factor(x);f[#f,1] отрабатывает сильно быстрее (менее сотой секунды). Без циклов и лишнего кода.

[wrest]

Я вот тоже что-то думаю, что быстрее всего будет как-то так, напрямую:
maxprimefactor(n)=my(f=factor(n));return(f[#f~,1])

-- 18.01.2024, 03:18 --

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

Но соотношения почти не изменятся.

Сдается мне, что поиск максимального простого делителя и полная факторизация -- одного поля ягоды (одной цены по сложности в среднем).
В отличие от поиска минимального (где перебор простых по возрастанию начиная с двойки будет в среднем быстрее).
P.S. Поэтому я весь свой первоначальный пост удалил, т.к. понял что фигню придумал

-- 18.01.2024, 03:28 --

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

if(x%p==0,print(x/1e15,"e15:",p)

Вам кстати тоже на заметку Если сравнивать [целое число] надо с нулём, то конструкция if(n==0,seq) работает заметно (на 10-15% по моим подсчетам) медленнее чем if(n,,seq)
Хотя конечно зависит от тяжести seq и частоты выпадения n<>0

[Dmitriy40]

Ну думаю да, практически факторизация, сразу об этом написал.
Про if видел, интересно.

[gris]

forvec я немного знаю, но нет ли там какого дополнения?
Дан вектор из векторов. Надо подать на обработку все вот такие векторы. У меня получается как-то запутанно. Может быть foreach приткнуть?
w=[[1,2],[3,4,5,6],[7,8,9]];
forvec(v=vector(#w,i,[1, #w[i]]),
print( vector(#w,j,w[j][v[j]])) )

[1, 3, 7],[1, 3, 8],[1, 3, 9],[1, 4, 7],
...
[2, 5, 9],[2, 6, 7],[2, 6, 8],[2, 6, 9]

[Dmitriy40]

Ну можно print(vector()) заменить на foreach(print()), но зачем, первое ведь короче и удобнее.

[gris]

спасибо. но вот вопрос: вот эти различные штучки красивы и удобны. Но не замедляют ли они вычисления? Ведь всё можно сделать через несколько for.
Если у меня намертво вбит массив w[345,554,787, ...] и его элементы многократно используются, то не эффективней ли везде вместо w[1] вручную написать 345 и т.д.
Вместо for( i=1,#d написать ld=#d; for( i=1,ld ? Если внешний цикл проигрывается триллион раз?
Я пробовал с таймером и в разных обстоятельствах, но неубедительно .

[wrest]

gris в сообщении #1626811 писал(а):

Я пробовал с таймером и в разных обстоятельствах, но неубедительно .

Ну тут теоретизировать дело неблагодарное, имхо, и лучше проверить. Но так-то да, статические константы теоретически работают быстрее, чем переменные. Компиляторы-интерпретаторы порой творят чудеса догадливости, а порой тупят, как дети.
Лично я всегда стараюсь, если возможно, сразу отделить длину (#) в переменную и использовать её. Но я не знаю точно что будет дольше, a=#v+1;b=#v*3;c=#v-1 или len=#v;a=len+1;b=len*3;c=len-1

[Dmitriy40]

gris в сообщении #1626811 писал(а):

не эффективней ли везде вместо w[1] вручную написать 345 и т.д.

Для PARI - надо проверять. Иначе никак, слишком он неоднозначный в плане скорости.
Но тут важнее другое: константы неудобнее человеку. Сильно неудобнее. Вдруг понадобится поменять длину массива или паттерн или ещё что? И ищи все эти числа по тексту ... И ещё не перепутай 1 как номер/индекс и 1 как величину числа и 1 как разницу между двумя объектами и 1 как код символа и 1 как булевкую истину и 1 как ещё что-нибудь ... Совсем не зря придуманы куча техник облегчения труда программисту, даже в ущерб скорости работы программы.
Хотите триллионы+ итераций - подумайте насчёт выучить другой язык, компилируемый, (или технологию gp2c для PARI), хоть С++, хоть Паскаль/Дельфи, хоть Фортран, хоть даже Бейсик - и переписать внутренние циклы на нём, это наверняка выйдет быстрее чем PARI досчитает.

wrest в сообщении #1626812 писал(а):

Лично я всегда стараюсь, если возможно, сразу отделить длину (#) в переменную и использовать её.

Неоднократно сталкивался когда вынос вычисления из цикла в переменную замедляет программу. А уж длину вектора и вовсе получить тривиально для PARI - она явно есть в объекте "вектор", интерпретация же что #x, что ld вряд ли различается по скорости (хотя лучше проверить).

-- 22.01.2024, 18:25 --

wrest в сообщении #1626812 писал(а):

Но я не знаю точно что будет дольше, a=#v+1;b=#v*3;c=#v-1 или len=#v;a=len+1;b=len*3;c=len-1

? for(i=1,10^7, a=#v+1;b=#v*3;c=#v-1; );
time = 11,279 ms.
? for(i=1,10^7, len=#v;a=len+1;b=len*3;c=len-1; );
time = 14,508 ms.

Вот так вот. Я ж говорю, с PARI многие стандартные техники оптимизации не работают.

[wrest]

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

Вот так вот. Я ж говорю, с PARI многие стандартные техники оптимизации не работают.

Ну да... Хотя если вынести len=#v за цикл, то будет незначительно быстрее. Но настолько, что оно того не стоит в плане оптимизации.

[gris]

Повадился упаковывать циклы в вектор. Но вот не уразумею такого. Допустим надо от простого числа отступить по 10 шагов вправо-влево.
{p=71;
v=vector(19);
v[10]=p;
for( i=1,9,
v[10+i]=nextprime(v[ 9+i]+1);
v[10-i]=precprime(v[11-i]-1);
);print(v)}
[31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109]
Как это сделать изящно, но без особого напряга.
+++ Dmitriy40, была такая идея, но я подумал, что вычисление номера простого числа это уж совсем большой напряг . Как и concat и сортировка двух векторов. Это и не изящно ничуть

[Dmitriy40]

Страшно медленно для больших чисел:
? p=71; k=primepi(p); primes([prime(k-9),prime(k+9)])
%1 = [31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109]

[gris]

Вот такая штучка. Для натурального числа n в десятичной записи вывести множество цифр в к-ичной записи с учётом возможных ведущих нулей, если длина записи меньше l.Типа такого
{k=6; l=7; n=49;
d=digits(n,k); S=Set(d); if( #d<l, S=setunion([0],S)); printf("%d_%d = %d_%d: %d\n", n, 10,d,k,S)}

49_10 = [1,2,1]_6: [0,1,2]
1343734_10 = [4,4,4,4,4,5,5,4]_6: [4,5]
5788765_10 = [5,8,5,4,5,13]_16: [0,4,5,8,13]

Нет ли чего красившее?

[wrest]

gris в сообщении #1628129 писал(а):

Нет ли чего красившее?

Чтобы вывод был красивше или скрипт?
Если цифр меньше 33, то можно же буквы использовать и слитно писать тогда а не массивом.
Вы напишите примеры, какой вы хотите "красивый" вывод...

У вас там кажется неправильно получается насчет ведущих нулей или я не понял задание.

gris в сообщении #1628129 писал(а):

5788765_10 = [5,8,5,4,5,13]_16: [0,4,5,8,13]