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

[wrest]

waxtep
Возможно, сначала преобразование вектора во множество (убирает повторения и сортирует вектор), а затем уже операции над множествами.
Два множества по 10000 элементов вычлись мгновенно
? sv1=Set(vector(10000,i,random()));
? ##
*** last result computed in 4 ms.
? sv2=Set(vector(10000,i,random()));
? ##
*** last result computed in 5 ms.
? sv3=setminus(sv1,sv2);
? ##
*** last result computed in 0 ms.
?

[Dmitriy40]

Посортируйте вектора (vecsort), а потом setminus отрабатывает за доли секунды даже на пятимиллионных векторах:

? default(parisize,6*10^8)
  ***   Warning: new stack size = 600000000 (572.205 Mbytes).
? a=primes([1,10^8]); b=vector(#a,i,i*6+1);
time = 922 ms.
? c=setminus(b,a);
time = 281 ms.
? #a
%1 = 5761455
? #c
%2 = 4700720
?

[waxtep]

Ой-ой-ой, отставить панику, это я один из векторов (совершенно неожиданно) формирую за какое-то безумное время. Сам расчет элементов менее минуты, но созданный мною наворот из concat и if почему-то безумно замедляет. Пока вопрос снимаю

-- 30.03.2022, 18:57 --

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

отрабатывает за доли секунды даже на пятимиллионных векторах

Да, я что-то не то делаю при формировании одного из векторов, пардон, буду разбираться.

-- 30.03.2022, 19:15 --

Вот в чем закавыка: закомментированный вариант функции sm3 безумнейше тормозит:

Используется синтаксис Matlab M

sm(n)={if(n<2, return(0)); while(1-isprime(n), n=n-numdiv(n)+1); return(n)};

\\sm3(s,f)={ng=matrix(1,2); for(k=s, f, if((g=sm(k))<k, ng=concat(ng,[k,g]))); return(ng)};

sm3(s,f)={ng=matrix(f-s+1,2); for(k=s, f, ng[k-s+1,1]=k;); for(k=s, f, if((g=sm(k))<k,ng[k-s+1,2]=g;)); return(ng)};

sm4(s,f)={ng=sm3(s,f); pp=Set(ng[,2]); return(pp)};

sm5(s,f)={p0=primes([1,f]); pp=sm4(s,f); sgp=setminus(p0,pp) ; return(sgp)};

В то время как раскоменнтированый, хоть и выглядит несимпатично, отрабатывает быстро и дальше проблем тоже нет. В чем здесь некорректность, работа с неинициализированной памятью, или какая-то более банальная ошибка?

-- 30.03.2022, 19:19 --

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

Два множества по 10000 элементов вычлись мгновенно

Да, не в том месте проблема, а ранее, тоже честно говоря для меня загадочная (дополнил свой последний комментарий)

-- 30.03.2022, 19:45 --

Или, если без попробовать без индусского кода, как-то так:

Используется синтаксис Matlab M

sm(n)={if(n<2, return(0)); while(1-isprime(n), n=n-numdiv(n)+1); return(n)};

sm5(s,f)={p0=primes([1,f]); pp=Set(vector(f-s+1,k,if((g=sm(s+k-1))<s+k-1,g,0))); sgp=setminus(p0,pp) ; return(sgp)};

[waxtep]

Возник такой еще вопрос, пытаюсь запустить один из примеров рисовалки из хэлпа, но получаю ошибку:

(01:35) gp > ploth(X=-1,1, sin(1/X), "Recursive")
  ***   at top-level: ploth(X=-1,1,sin(1/X),"Recursive")
  ***                              ^---------------------
  *** sin: precision too low in mpcosm1.

Насколько я смог понять, mpcosm1 - внутренняя функция, на которую опирается вычисление косинусов и синусов, но... кхм как с ней бороться и что я делаю не так, - непонятно

-- 31.03.2022, 01:52 --

Думал, что проблема в окрестностях нуля, но нет - около единицы почему-то:

(01:50) gp > ploth(X=-1,1-1e-19, sin(1/X), "Recursive")
%16 = [-1.0000000000000000000, 1.0000000000000000000, -0.9999985392207986212, 0.9999985392207986212]
(01:50) gp > ploth(X=-1,1-1e-20, sin(1/X), "Recursive")
  ***   at top-level: ploth(X=-1,1-1e-20,sin(1/X),"Recursive")
  ***                                    ^---------------------
  *** sin: precision too low in mpcosm1.

[Dmitriy40]

waxtep
Вызываете синус со слишком большим аргументом, и PARI не хватает точности представления чисел чтобы получить остаток от деления на $2\pi$.
И это не лечится.
Исключите слишком близкие к нулю координаты по X (ближе чем примерно $10^{-30}$ для стандартной точности в 38 цифр). Например так:

? ploth(X=-1,1, if(abs(X)>10^-30, sin(1/X), 0.1), "Recursive")
[-1.0000000000000000000, 1.0000000000000000000, -0.9999985392207986212, 0.9999985392207986212]

Или перейдите в поле рациональных чисел (не знаю как), в нём синус должен вычисляться точнее.

Проблема около нуля: когда оно остановилось по ошибке наберите X и посмотрите чему он равен:

? ploth(X=-1,1-10^-29, sin(1/X), "Recursive")
  ***   at top-level: ploth(X=-1,1-10^-29,sin(1/X),"Recursive")
  ***                                     ^---------------------
  *** sin: precision too low in mpcosm1.
  ***   Break loop: type 'break' to go back to GP prompt
break> X
-5.000000001161704870 E-30
break> 1/X
-1.9999999995353180522 E29
break>

[waxtep]

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

Проблема около нуля: когда оно остановилось по ошибке наберите X и посмотрите чему он равен:

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

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

*** Break loop: type 'break' to go back to GP prompt
break> X

А как Вы такую штуку делаете? Это надо какой-то параметр отладки выставить?

-- 31.03.2022, 02:12 --

ploth(X=-1,1, if(abs(X)>10^-30, sin(1/X), 0.1), "Recursive")

О, спасибо! Так работает

[Dmitriy40]

waxtep в сообщении #1551473 писал(а):

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

Вы ошибаетесь, беда около нуля, за ноль в положительные числа он вообще не заходит. А то что Вы туда добавили вещественное число лишь чуточку сместило точку возникновения ошибки.

waxtep в сообщении #1551473 писал(а):

А как Вы такую штуку делаете? Это надо какой-то параметр отладки выставить?

У меня изначально так работает. Но вообще да, выставьте default(breakloop,1) и будет Вам счастье (кстати сам не знал про эту возможность).
Текущие параметры можно посмотреть по \d.

-- 31.03.2022, 02:21 --

waxtep в сообщении #1551473 писал(а):

О, спасибо! Так работает

Зацените график если вместо 0.1 поставить просто 0 ... О причинах не спрашивайте, не знаю.

[waxtep]

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

выставьте default(breakloop,1) и будет Вам счастье

Йес! Спасибо

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

Зацените график если вместо 0.1 поставить просто 0 ... О причинах не спрашивайте, не знаю.

Жуть ну, они пишут, что "Although plotting is not even a side purpose of PARI..."
Кстати, эту беду можно еще вот так победить, поменять количество точек для рекурсивного рисунка на нечетное, тогда мы видимо успешно промахиваемся мимо нуля:

(02:25) gp > ploth(X=-1,1, sin(1/X), "Recursive",7)
%1 = [-1.0000000000000000000, 1.0000000000000000000, -0.9999955771745117072, 0.9999548190565817896]
(02:26) gp > ploth(X=-1,1, sin(1/X), "Recursive",8)
  ***   at top-level: ploth(X=-1,1,sin(1/X),"Recursive",8)
  ***                              ^-----------------------
  *** sin: precision too low in mpcosm1.

-- 31.03.2022, 02:36 --

Вот она, попалась:

break> ploth(X=-1,1, sin(1/X), "Recursive",9)
  ***   at top-level: ploth(X=-1,1,sin(1/X),"Recursive",10)
  ***                              ^------------------------
  ***   in function sin: ploth(X=-1,1,sin(1/X),"Recursive",9)
  ***                                      ^------------------
  *** _/_: impossible inverse in divrr: 0.E-38.
  ***   Break loop: type 'break' to go back to GP prompt
break[2]> X
0.E-38

-- 31.03.2022, 02:44 --

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

Зацените график если вместо 0.1 поставить просто 0 ...

Это видимо вот что:

мануал писал(а):

For instance you should never plot recursively an odd function in a symmetric interval around 0. Try

ploth(x = -20, 20, sin(x), "Recursive")

[Dmitriy40]

waxtep в сообщении #1551475 писал(а):

break[2]>

Не забывайте теперь после каждой ошибки выходить на самый верхний уровень командами break, вот тут этого не сделали и запустили новую команду в контексте запущенной предыдущей. Могут появиться самые разные неприятные побочные эффекты.
А просто нажатие Enter приведёт к продолжению выполнения прерванной (например по Ctrl+C) программы/функции. Например можно поменять какие-то переменные и продолжить выполнение дальше.

[waxtep]

Dmitriy40, спасибо! Будем ломать и строить

[Dmitriy40]

Камрады, кто-то может пояснить что за хрень иногда происходит и как её побороть?

  ***   at top-level: .................;f=alarm(1,factor(x));
  ***                                             ^----------
  *** factor: user interrupt after 1h, 1min, 57,270 ms
  ***   Break loop: <Return> to continue; 'break' to go back to GP prompt
break>
  ***   at top-level: ..............;alarm(3,s[d]=numdiv(x));
  ***                                             ^----------
  *** numdiv: user interrupt after 3h, 35min, 58,240 ms
  ***   Break loop: <Return> to continue; 'break' to go back to GP prompt
break>

Если текст вдруг съедет, то стрелка указывает на factor и numdiv.

В упор не понимаю почему user interrupt, если никакой юзер там и близко ничего не делал, а вылетает программа сама и на разных числах. Причём по Enter вполне нормально продолжает считать дальше, т.е. ошибка не критическая. И по всему ошибка должна быть e_ALARM, а вовсе не e_USER (e_ALARM я дальше сам обрабатываю).
Если в остановившимся PARI руками ввести ровно ту же команду factor или numdiv вместе с alarm, то возникает совершенно правильная e_ALARM и никаких user interrupt нет и в помине.
Возникает спонтанно, повторный запуск тех же вычислений проходит нормально.
От версии PARI не зависит, проверил и на самой последней.

[EUgeneUS]

Dmitriy40
Что у Вас с настройками сна и гибернации? Не связан ли "user interrupt" с попытками компа уйти в сон или гибернацию?

(Оффтоп)

при запуске скриптов от VAL столкнулся с подобным:
1. При переходе в сон комп действительно засыпал. Но при выходе из сна - счет продолжался.
2. (Предположительно) при переходе в гибернацию счет ВНЕЗАПНО завершался. Так как PARI/GP запускал через ярлык, окна с командной строкой закрывались, и диагностика терялась.
Что интересно, при расчете 15-12 с использованием "ускорителей" такого не было (комп не пытался уснуть или впасть в гибернацию, хотя настройки такие были). Предполагаю это связано с тем, что в этом случае происходило регулярное обращение к HDD.

[Dmitriy40]

Ни в сон ни в спячку сам не уходит, он работает 24/7/365 (кроме вычислений крутятся ещё другие сервисы).
При останове в одном потоке в остальных счёт продолжается (пока тоже не вылетят, например за ночь вылетели три потока из 4-х, с разницей более чем в 5ч).
Главное не понимаю почему не ту ошибку вываливает, или даже вообще не ошибку (ошибки я как раз перехватываю и обрабатываю). Смахивает на глюк CPU/RAM или самого PARI.

Бесит меня этот глюк, один поток 11ч проработал, а другой трижды за час вылетел.

По скриптам VAL — так там же включен лог, в нём всё должно быть.

[dmd]

Dmitriy40
Может включение дэбага подскажет что происходит: ключи \g, \gm, \gf, \h?

[Dmitriy40]

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