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

[Dmitriy40]

Тут опять всё зависит от длины числа в битах и от количества единичных битов. Если единичек несколько штук, то Ваш вариант возможно будет быстрее, если же единичек [может быть] много, то я бы сделал так:

? v=binary(n); len=matsize(v)[2]; for(i=1, len, if(v[i]==1, <что делать при единичках>));

Перебор всех 10 тысяч единиц в числе $2^{10000}-1$ занимает у меня около 70мс. А сто тысяч единичек в числе $2^{100000}-1$ за 270мс.

Да, если так делать, не забудьте что биты тут нумеруются параметром i со старшего и с 1, а не с 0.

[dmd]

Кстати, получить вектор индексов установленных битов можно так:

select(i -> i==1, binary(n), 1)

-- Чт июн 02, 2016 00:28:36 --

Вопрос. Как написать свою функцию (или процедуру) forbit(n, X, sec), которая была бы подобна встроенной fordiv(n, X, sec)?

[dmd]

Снова вопрос про решение уравнений по составному модулю. Допустим уравнение кубическое по модулю 9.


Mathematica показывает три правильных решения:

Reduce[2 + X*4*(3 + X*4*(3 + X*4)) == 0, {}, Modulus -> 9]
X == 1 || X == 4 || X == 7

pari/gp по простому модулю 3 показывает одно:

polrootsmod(2+X*4*(3+X*4*(3+X*4)), 3)
[Mod(1, 3)]~

Как в pari/gp получить все три решения по модулю 9?

[maxal]

dmd в сообщении #1157351 писал(а):

Как в pari/gp получить все три решения по модулю 9?

Разве что "подъёмом" по лемме Гензеля:

? ?polhensellift
polhensellift(A, B, p, e): lift the factorization B of A modulo p to a factorization modulo p^e using Hensel lift. The factors in B must be pairwise relatively prime modulo p.

[Dmitriy40]

А получится ли? Если для $\mod 9$ PARI ещё выдаёт корень $(1)$, к тому же все три корня равны по $\mod 3$, то для $\mod 27$ PARI выдаёт пустое множество корней, а WolframAlpha - три корня $(4, 13, 22)$. Получится ли лифтом из пустого множества получить корни? К сожалению не понимаю как самому проверить.

-- 05.10.2016, 18:32 --

Для $\mod 81$ и $\mod 243$ ситуация аналогична.
Кстати, в polrootsmod есть третий параметр, если его указать в 1, то для модулей 27, 81 и 243 появляется корень, но снова лишь один наименьший.

[maxal]

Dmitriy40, polrootsmod определён только для простых модулей. Соответственно, на непростых модулях его результат неопределен, и лучше на него не полагаться.

-- Wed Oct 05, 2016 10:54:05 --

P.S. Лемма Гензеля работает только для свободных от квадратов многочленов (где производная на корнях не обращается в нуль) по простому молулю. Про сложности с кратными корнями см. в википедии.

[Dmitriy40]

Ясно, спасибо, не подумал что простота принципиальна, всё, больше не мешаю.

[Коровьев]

Мне эта программа очень нравится. К сожалению, программирование моё застыло на началах Бейсика, и по сему много вопросов по PARI/GP.
Вот вопрос.
При длительных вычислениях абсолютно не ясно, то ли идти спать, или ещё потаращиться на экран в ожидании конечного результата.
В Бейсике, помню, я расставлял "вехи" при переборе переменных, к примеру -
идет идёт последовательный перебор двух переменных $x,y$. Когда меняется вторая переменная $y$ следует команда "PRINT $y$. Но каждый раз печать идёт на одно и то же место экрана. В Бейсике есть такая команда, поэтому всё время на экране меняется только одна строчка. По ней я судил где находится программа вычисления и надо ли идти спать..
Есть ли такая команда в PARI/GP?
Ест ли другие возможности контроля хода программы вычислений?

[maxal]

Коровьев
Вот примерный аналог:

print1(y," ",Strchr(13));

[Коровьев]

maxal
Большое спасибо! Всё работает.
Вот ещё у меня проблема.
Определить является ли квадратом полученное целое число довольно просто.

k=abs(x)-(sqrtint(abs(x)))^2

Если $k=0$, то $x$ есть квадрат целого числа.
Но как командно определять: являются ли квадратом рационального числа получаемые неким перебором дробные числа? Приходится изворачиваться.
И вот вопрос.
Есть ли возможность командно из дроби выделять числитель (или знаменатель)?

[maxal]

Коровьев в сообщении #1175579 писал(а):

Определить является ли квадратом полученное целое число довольно просто.

Для этого есть функция issquare(). Работает как для целых, так и для рациональных чисел.

Коровьев в сообщении #1175579 писал(а):

Есть ли возможность командно из дроби выделять числитель (или знаменатель)?

numerator(), denominator()

[Aether]

Вот код, который находит последовательность чисел, являющихся целыми степенями чисел Мерсенна:

ismn(n) = n++; n == 2^valuation(n, 2);
isok(n) = ismn(n) || (ispower(n, , &m) && ismn(m));

Код очень, просто удивительно короткий.

Нельзя ли подробно разобрать этот код?

[arseniiv]

(Комментарий от не знающего PARI/GP)

По поводу valuation см. https://math.mit.edu/~brubaker/PARI/PARIrefcard.pdf (что она делает) и https://en.wikipedia.org/wiki/P-adic_order (что это значит), так что ismn определяет, является ли её аргумент числом Мерсенна, и вполне прямым образом.

Вот по поводу второй функции, а особенно выражения ispower(n, , &m), уже нужно разбираться в PARI/GP.

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

который находит последовательность чисел

Да ну, просто же проверяет, такое число или нет.

[Aether]

Хорошая и познавательная табличка по PARI\GP, большое спасибо. Думаю, её неплохо было бы выложить вначале темы.

В OEIS код заявлен как строящий последовательность.

[arseniiv]

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

Хорошая и познавательная табличка по PARI\GP большое спасибо? Думаю её неплохо было бы выложить вначале темы.

Она 2003 года для версии 2.2.5. Думаю, там с тех времён уже много изменилось. В частности, я в ней не нашёл ispower. Также, maxal в первом посте темы использовал 2.4.3, а сейчас доступна уже 2.9.1.

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

В OEIS код заявлен как строящий последовательность.

Ну, если множество разрешимо (есть вычислимая функция, определяющая, принадлежит ему элемент или нет), то оно и перечислимо (есть вычислимая функция, «перечисляющая» все его элементы — реализующая произвольную его нумерацию), потому что мы можем перебирать в цикле все натуральные числа и выдавать те из них, которые удовлетворяют предикату. Написать такую «обёртку» не представляется сложным на типичном языке программирования.