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

[Dmitriy40]

За 5ч уже полгига лога, а глюка так и нет.
Причём в лог явно попадает не всё, например никак не отражаются вызовы ispseudoprime, хотя поставил максимальный уровень дебага. Они конечно и не интересны, но непонятка.
Что полезного смотреть в \h вообще не понял.

[Dmitriy40]

А вот и глюк.
Но в логе полная тишина о причинах, ровно то же уже процитированное сообщение:

...
ECM: finishing curves 0...3
        (extracted precomputed helix / baby step entries)
        (baby step table complete)
        (giant step at p = 28697)
        (giant step at p = 82457)
        (giant step at p = 136217)
ECM: time =    110 ms
ECM: B1 = 2000, B2 = 220000,    gss =  128*420
  ***   at top-level: .................;f=alarm(1,factor(x));
  ***                                             ^----------
  *** factor: user interrupt after 4h, 32min, 46,002 ms
  ***   Break loop: <Return> to continue; 'break' to go back to GP prompt
break>

Время показано меньше реального, это особенность PARI с подсчётом времён.
Те строки что прямо перед ним встречались за эти 6ч тысячи раз, что-то там про разложение на простые.
Всего конструкция alarm(1,factor()) до появления глюка сработала правильно примерно 8000 раз.

Выходит это натуральный глюк в самом PARI в работе команды alarm, потому что без alarm ни factor ни numdiv до сих пор таких глюков не давали.

[dmd]

Все-таки неисправность железа тоже равновероятна.

(Оффтоп)

Быть может в одном из физических ядер проца произошла деградация, и желательно его отключение (если биос конечно позволяет). Как тестировать ядра вообще не представляю. Либо память. Протестировать ее легко, но муторно долго. Утилитой TestMem5 снaчала с конфигом "default" на два-три прохода (часа 2 где-то длится). Затем с конфигом "[e-mail]" на один-два прохода (сутки-двое будет по времени, зато точно станет понятно выдерживает ли память такую нагрузку).

[Dmitriy40]

Ядра можно проверить так: назначить каждому потоку своё ядро и смотреть где будут вылеты. Пожалуй сделаю тестовую программку и проверю, чисто ради интереса.
Память проверял недавно, пару месяцев назад, тоже было подозрению на странные глюки, нет, всё в порядке.
Но сильно сомневаюсь что дело в аппаратуре: тогда глюки были бы и с другими вычислениями (в том числе и ровно таких же factor/numdiv), а их не было и нет. И нет не потому что не поймал, как-то раз комп всю ночь и полдня считал уже насчитанное другим компом и в итоге всё совпало. Вот и сейчас за ночь (13ч счёта) 4 потока ни разу не сбились, хотя всего лишь убрал alarm, остальной код фактически тот же.

[aitap]

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

*** at top-level: .................;f=alarm(1,factor(x));
***                                           ^----------
*** factor: user interrupt after 4h, 32min, 46,002 ms

PARI работает на Windows в однопоточном режиме? Думаю, единственный способ получить такую ситуацию - это попасть в gp_sigint_fun (которая, судя по всему, используется на Windows для обработки alarm(s, code)), но таким образом, что переменная win32alrm почему-то не выставлена.

Возможно ли это? win32alrm выставляют в win32_cb_alarm, причём сначала выставляют флаг win32ctrlc (должно быть, значит "нужно запустить gp_sigint_fun"), и только потом включают win32alrm. Как запускается эта функция? Чуть ниже, её готовят для запуска при помощи CreateTimerQueueTimer с флагом 0. By default, the callback function is queued to a non-I/O worker thread.

Думаю, этот рабочий поток успевает выставить win32ctrlc, после чего управление достаётся PARI, которая его проверяет и запускает обработчик user interrupt, а win32alrm остаётся невыставленным на момент его проверки основным потоком. Либо потоки работают на разных ядрах, и в отсутствие какой-либо синхронизации между ними записи одного потока в память не сразу становятся видны из другого потока. Я даже не скажу навскидку, как это исправить; по-честному, тут нужны мьютексы или atomics, чтобы убедить и компилятор, и процессор не менять записи в память местами в целях оптимизации. Ну и, по меньшей мере, строчку с записью в win32alrm нужно поставить до записи в win32ctrlc.

[Dmitriy40]

aitap в сообщении #1553229 писал(а):

PARI работает на Windows в однопоточном режиме?

Да. Многопоточность в PARI очень уж кривая (вместо скорости х4 стабильно выдаёт лишь х3 при 100% загрузке всех 4-х ядер).
Спасибо за разбор! Давненько не лазил в исходники PARI. Не очень понял возможна ли эта ситуация в однопоточном режиме (ведь callback функция будет вызвана не асинхронно в момент окончания времени таймера, а синхронно при следующем опросе очереди сообщений, как мне представляется). Ну и выходит без перекомпиляции самого PARI это не обойти ...

[wrest]

(Оффтоп)

Dmitriy40
Истинно говорю вам: да перейдъте уже как-нибудь на WSL или Линукс с minGW...

[aitap]

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

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

Думаю, возможна. Судя по документации CreateTimerQueueTimer, Callback запускается в отдельном потоке. ("Callback functions are queued to the thread pool.")

[Dmitriy40]

Кажется наткнулся на решение этой проблемы, во всяком случае пока вылетов больше не было (правда и проверок было на пару порядков меньше): достаточно перед alarm(1,...) сделать alarm(0). Наверное можно и после, но менее удобно, там же ветвление.
В принципе это даже в мануале приведено, только слишком туманно разъяснён смысл (сбросить флаги сработавших таймеров). Но если и правда работает, то выходит глюк возникал из-за срабатывания таймера на предыдущий вызов alarm(), который не отменялся ни при выходе из предыдущего, ни при входе в новый, а сразу срабатывал. Или сначала проверялся, а уж потом перезапускался. И пока между вызовами alarm(1,...) проходило менее секунды всё работало, а как достигало секунды, то следующий alarm() вылетал.
Хотя нет, тест эту версию не подтверждает:

x=16774592909948271587441923785827599104180389597956520637078459491512670623892264042491;
alarm(1,print(numdiv(x)));
n=0; forprime(p=1,10^9, n++);
alarm(1,print(numdiv(x)));

Вполне работает, хотя цикл занимает сильно больше секунды.
Значит возможно alarm(0) проблему и не лечит, просто ещё не наткнулся на глюк.

[Dmitriy40]

Да, alarm(0) не лечит, вот только что глюк возник и с ним.

[dmd]

Возможно ли получить сумму делителей быстрее чем sumdiv(n) в однопотоке и D=divisors(n); export(D); parsum(j=1, #D, D[j]) в многопотоке? Вопрос возник в связи с задачкой, для которой уже $n_{31}$ любое хотя бы не рекордное у меня не получается найти.

[Dmitriy40]

dmd
У Вас же не произвольное число чтобы использовать sumdiv, замените его на вычисление суммы делителей прямо по вектору простых:

? pr=setminus(primes([17,50]),[19,37,41]); print(sumdiv(vecprod(pr),i,i)); print(prod(i=1,#pr, pr[i]+1));
875888640
875888640

А параллелить лучше не внутренность цикла (sumdiv), а сам цикл перебора.

[Dmitriy40]

dmd

(Оффтоп)

Собственно в той задаче и не нужна sumdiv, там перебор можно организовать по сумме векторов степеней в факторизациях. Вот к примеру найдено руками (по насчитанным таблицам) за пару часов:
p=31: [31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 109, 127, 137, 149, 163, 173, 181, 191, 197, 199, 257, 269, 311, 401]
p=37: [37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 109, 127, 137, 163, 173, 181, 191, 199, 233, 239, 257, 269, 271, 311, 367, 401, 449, 487, 521, 593, 599]

[dmd]

Dmitriy40

(Оффтоп)

Здорово! Дайте там свой коммент или ответ.
Сколько еще неожидаемого даст эта задачка? Меня впечатлило "квантование" при удалении любого простого из $n$, когда "след" удаленного остается виден в $k(n/p)$.

[Dmitriy40]

dmd

(Оффтоп)

Да там оказывается они и сами нашли и p=31 и p=37, мои варианты больше. Разве что найти p=41 и p=43 и p=47 ...
Зато кажется у меня немного другой метод нахождения, более общий, позволяет относительно быстро найти почти любое p, по заранее посчитанным общим таблицам (которые правда считать как раз долго), но решение будет далеко не минимальным.