08.10.2024 09:01
inetstar 18 1900 Источник


Мир изменился. Многоядерные процессоры повсюду. Использование их потенциала позволяет количество вычислительной мощи превратить в новое качество. Многие задачи стало выгодно численно решать и проводить эксперименты над ними.

В этой статье я расскажу о программе PARI, язык которой GP имеет очень низкий порог вхождения, примерно как Basic, но при этом позволяет проводить быстрые сверхточные целочисленные и вещественные вычисления (см. предыдущую статью). А большой свечкой на торте является то, что с её помощью можно совершенно элементарно загрузить все ядра ваших процессоров и серверов для получения максимально быстрого результата.

В этой статье за ~15 минут вы узнаете, как легко и просто загрузить компьютер на 100% вашими вычислительными задачами, даже если вы не являетесь профессиональным программистом.

Оглавление


0. Для кого эта статья
1. Что это за программа PARI/GP
2. Как писать программы для PARI
3. Особенности синтаксиса
  3.1 Локальные и глобальные переменные
  3.2 Блоки кода
  3.3 Каждый оператор может что-то возвращать
  3.4 Условный оператор
  3.5 Циклы
  3.6 Создаём массивы, списки, множества и словари
  3.7 Передача по ссылке и указателю
4. Параллельные вычисления
  4.1 Особенности параллельных вычислений
  4.2 Конструкция parfor и другие
  4.3 Уникальные случайные числа в каждом треде
  4.4 Компиляция медленных функций
  4.5 Другие параллельные конструкции
5. Хитрости, которые я использую в Linux
  5.1 Утилита sort
  5.2 Утилита uniq
6. Итоги

0. Для кого эта статья


Вам нужно быстро посчитать что-то математическое/криптографическое (тут пригодятся масса математических функций и целые неограниченного размера) или физическое (тут пригодятся
float неограниченной точности), желательно параллельно (для ускорения работы) и с минимальным порогом вхождения (через полчаса вы уже получите работающий скелет параллельной программы или её всю), то вам сюда. Будем изучать программу PARI/GP, которую разработчики когда-то назвали «Great Programmable Calculator», но она выросла в нечто значительно бОльшее.

1. Что это за программа PARI/GP


PARI/GP — это система компьютерной алгебры, в первую очередь предназначенная для специалистов по теории чисел, но из-за своей скорости, точности и простоты нашедшая применение во многих других областях, от топологии, криптографии или численного анализа до физики. HTML Docs, Tutorial, User guide, шпаргалка.

PARI умеет делать символьные манипуляции, но справляется с такими задачами хуже, чем специализированные CAS по типу Maxima, Maple или Mathematica. С другой стороны, три основных преимущества системы — это её скорость, множество специальных типов данных, которые знакомы математикам, и богатое количество функций из теории чисел.

Я пишу статью о PARI/GP из-за того, что был поражён простотой входа в сочетании с вычислительными возможностями программы, в том числе параллельными.

Нематематические сильные стороны:

  1. Высокоуровневый скриптовый язык.
  2. Программирование на языках общего назначения с помощью библиотеки PARI.
  3. Надёжность и зрелость (разработка началась в 1985 г.).
  4. Хорошая поддержка для пользователей со стороны разработчиков (например, ответили на моё письмо через несколько часов).

И, конечно же, PARI/GP является свободным программным обеспечением, на которое распространяется лицензии GNU General Public License v.2 или более поздняя.

PARI используется тремя различными способами:

  1. Как сложный программируемый калькулятор, называемый gp, язык GP которого состоит из типовых управляющих конструкций, похожих на язык C.
  2. Как библиотека libpari, которую можно вызвать из приложения на языке верхнего уровня, например, написанного на C или C++.
  3. Компилятор gp2c транслирует код GP в C и загружает его в интерпретатор gp. Типичный скрипт, скомпилированный gp2c, работает в 3–10 раз быстрее. Сгенерированный код C можно редактировать и оптимизировать вручную.

Самый простой способ попробовать PARI/GP — это браузерная версия (скомпилированный WebAssembly).

Скачать версии под Linux, MacOS и Windows.

2. Как писать программы для PARI


Программы можно писать в интерактивном режиме, если нужно что-то совсем простое посчитать.

Вот как выглядит gp, запущенная в интерактивном режиме:



Либо в виде скриптов *.gp. Для запуска скрипта я использую команду:

gp -q script.gp

Если во время исполнения возникает ошибка, то gp переходит в отладочный режим и вы можете вывести состояние переменных, что-то исправить и продолжить выполнение.

Для того, чтобы выйти из программы, можно использовать quit() или \q.

Ещё полезные команды:

  • ## — вывод времени, потраченного на последнюю команду;
  • ?? — открытие в отдельном окне pdf-файла документации.

gp при запуске резервирует себе область памяти, в которой работает так называемый стэк. По умолчанию он около 800КБ и его может не хватить. Для его увеличения можно использовать команды:

\\ удвоение стека
allocatemem()
  ***   Warning: new stack size = 16000000 (15.259 Mbytes).
\\ Или задаём конкретный размер стека
default(parisize, 10*10^6);
  ***   Warning: new stack size = 10000000 (9.537 Mbytes).

Как вы заметили, для однострочных комментариев используются два обратных слэша. А для многострочных — стандартная конструкция /*… */.

Каждое выражение может что-то возвращать, точка с запятой в конце выражения подавляет этот вывод. Это в интерактивном режиме.

? x=1
%1 = 1
? x=1;
?

Что приятно — PARI работает с естественным порядком математических операций, и общепринятыми приоритетами операндов.

? 4+3*2^3-1
%3 = 27


3. Особенности синтаксиса


▍ 3.1 Локальные и глобальные переменные


Для объявления локальных переменных функций используется конструкция my(). Объявленные переменные можно использовать для объявления следующих переменных. Пример:

test(x) =
{
  my(
    a = 1,
    b = a/3
  );
  return (x+a+b)
}

▍ 3.2 Блоки кода


Как было видно из предыдущего примера, для объявления функции мы использовали конструкцию из фигурных скобок. Но в принципе мы могли бы задать функцию в одну строку.

? text(x) = my(a=1, b=a/3); return(x+a+b);
? text(1)
%29 = 7/3

Видно, что по возможности PARI старается дать точный ответ. Тип этого ответа t_FRAC (дробь). Если нам нужен вещественный ответ, то нужно было бы, чтобы одно из чисел было бы с точкой. Или написать return(x+a+b+0.0).

? text(x) = my(a=1, b=a/3); return(x+a+b+0.0);
? text(1)
%32 = 2.3333333333333333333333333333333333333

Для склеивания операторов в один блок применяется ;. Также ; используется для разделения операторов, как в С.

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

А запятая используется для отделения аргументов. Многие операторы в PARI Lisp-подобные, они выглядят как функции. В этом нет ничего сложного.

Зачем вообще склеивать операторы?

Рассмотрим пример условного оператора:

if (a > 2,
  \\ Если условие истинно
  print("a > 2"), \\ Запятая показывает, что ветка истины закончилась

  \\ Если условие ложь
  print("a > 2");
  a *= 2
)

▍ 3.3 Каждый оператор может что-то возвращать


Мы помним, что операторы в PARI как функции. Поэтому не удивительно, что они что-то возвращают. Это можно использовать, а можно не использовать.

? b = if(3 > 2, 3, 2);
? b
%34 = 3

▍ 3.4 Условный оператор


Кроме условного оператора с двумя ветками, мы можем использовать и более простой только с веткой «истины». Он иногда более удобен, так как позволяет не громоздить внутри условного оператора сложную логику.

example_abs(x) =
{
  if(x > 0, return(1));
  if(x < 0, return(-1));
  0
}

Как видим, в конце мы обошлись без return(). Последний оператор возвращает ноль в вышестоящий блок кода.

▍ 3.5 Циклы


Циклов в PARI существует великое множество. Включая такие экзотические, как перебор по разным комбинациям, размещениям, простым числам, элементам массивов. Самый базовый и простой цикл выглядит так:

for(i=1, 10,
  print(i)
)

Переменная i существует только внутри цикла.

Для того, чтобы начать цикл с начала используется команда next() (аналог continue из других языков). Для выхода из цикла можно использовать break(число_циклов).

▍ 3.6 Создаём массивы, списки, множества и словари


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

Создадим массив с квадратами чисел от 1 до 5. Можно сделать так:

v = [1, 4, 9, 16, 25]

А можно так:

v = vector(5, i, i^2)

Подобных инициализаторов в PARI довольно много на любой вкус, что позволяет очень компактно объявлять и инициализировать массивы.

Можно объявлять массивы и с помощью интервалов:

? [-2..3]
%2 = [-2, -1, 0, 1, 2, 3]

vector() — это функция, возвращающая массив. Но есть ещё и конструкторы Vec(), Set(), List(), Map(). Для создания, соответственно, массивов (векторов), множеств, списков и словарей (ассоциативных массивов).

Из-за того, что у вектора (массива) не может быть изменён размер, существует тип список. В нём можно динамически удалять/добавлять элементы.

? l = List([4, 2, 1, 0]);
\\ Вставить в конец
? listput(~l, 7); l
%15 = List([4, 2, 1, 0, 7])
\\ Вставить на первое место, остальное сдвигается вправо
? listinsert(~l, 8, 1); l
%16 = List([8, 4, 2, 1, 0, 7])

Список также поддерживает индексную адресацию через квадратные скобки. Конструктор List() в качестве аргумента принимает вектор.

? l = List([-2..5])
%43 = List([-2, -1, 0, 1, 2, 3, 4, 5])

Интересно, что сортировка списков происходит немного быстрее, чем векторов. Это происходит потому, что списки сортируются по ссылке.

? l = List([4, 2, 1, 0]);
? listsort(~l);
? l
%4 = List([0, 1, 2, 4])

Сортировка массива:

? v = [4, 2, 1, 0];
? v1 = vecsort(v);
? v
%7 = [4, 2, 1, 0]
? v1
%8 = [0, 1, 2, 4]

Также поддерживаются интервалы в индексах.

vector(10, i, i)[1..3]
%45 = [1, 2, 3]

В типе множество каждый элемент может присутствовать только один раз.

? Set([1,1,2,2])
%44 = [1, 2]

Создаём словарь, меняем элемент, выводим элемент, проверяем существование значения по ключу:

?  M = Map(["a",1; "b",3; "c",7]); M
%4 = Map(["a", 1; "b", 3; "c", 7])
? mapput(M, "a", 2); M
%5 = Map(["a", 2; "b", 3; "c", 7])
? mapget(M, "a")
%6 = 2
? mapisdefined(M, "d")
%7 = 0


▍ 3.7 Передача по ссылке и указателю


Для передаче по ссылке тильду ~ нужно использовать дважды:
  1. в объявлении функции перед именем переменной
  2. при вызове функции перед именем переменной.

По ссылке можно передать только массив (вектор), список и матрицу (тип, который я не осветил в этой статье).

По указателю можно передавать любую переменную, но только в документированные функции PARI допускающие это, например в mapisdefined(M, x, {&z}). Самому в GP нельзя создать функцию с указателем. Фигурные скобки означают, что параметр необязателен. При передаче в функцию такого указателя на переменную в переменную будет что-то записано. При вызове функции также указание & обязательно.

4. Параллельные вычисления


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



Обратите внимание на последнюю строчку. То, что внутри компьютера выполнялось 3839 часов, заняло реального времени всего 65 часов.

И тут то же самое:


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

Многие задачи при такой вычислительной мощи стало проще решить численно или с помощью математического моделирования, или, например, вероятностными методами (Монте-Карло и аналогичными).

PARI может запросто загрузить все логические ядра вашего процессора. И так и сделает, если вы не уменьшите число доступных ему процессоров.

Однако если вы не используете специальных параллельных конструкций, то беспокоиться об этом не стоит. Так как немногие операторы GP поддерживают внутренний параллелизм.

Так мы разрешим использовать все логические процессоры, кроме двух.

default(nbthreads, default(nbthreads)-2);

Разработчики написали введение в параллельное программирование на PARI.

▍ 4.1 Особенности параллельных вычислений


Конечно, PARI не поддерживает такую гибкость, как Go в области параллельных вычислений, но для обычных математических/физических вычислений возможностей достаточно.

Ограничения такие, что функции, запущенные в параллельных тредах, не могут устанавливать глобальные переменные. Такие функции должны быть использованы как вещи в себе, без сайд-эффектов, как говорится.

Но читать глобальные переменные они могут, но не все, а только те, что передаются в параллельные процессы с помощью оператора export(var1, var2, var3).

Можно передать вообще все глобальные переменные в параллельные треды с помощью exportall(). Но я не рекомендую чисто с точки зрения контроля.

▍ 4.2 Конструкция parfor и другие


Язык GP поддерживает множество конструкций для параллельной работы, но я рассмотрю только одну, базовую, которой достаточно для 99% случаев.

Это параллельный цикл parfor(). Пример:

\\ Выведем все числа от 1 до 1000000, квадраты которых минус (число+1) есть простые числа.
testpar(n) =
{
  my(n=0);
  parfor(i=1, n,
    \\ выражение для параллельного исполнения
    i^2-(i+1),
    \\ результат выражения записывается в r
    r,
    \\ критическая секция, эта секция не может выполняться параллельно
    \\ результаты выстраиваются в очередь сюда на обработку
    if (isprime(r), print(i))
  )
  print("Num of ", n)
}

? testpar(10^6)
..
924209
858390
20
21
..
Num of 139483
\\ Мы что, случайно нашли способ сгенерировать числа с повышенным содержанием простых? Их около 14%!

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

Обычно в качестве выражения используется функция. Это удобно и тем, что функцию позднее можно скомпилировать в нативный код для ускорения.

Переменные i и r локальны в пределах parfor. И они передаются в критическую секцию, что позволяет связать результат и номер итерации.

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

▍ 4.3 Уникальные случайные числа в каждом треде


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

Я придумал для этого инициализировать начальный параметр генератора через /dev/urandom (linux-специфичный способ).

\\ Мы считываем целое беззнаковое число до 2^64 из /dev/urandom
setrand(extern("od -vAn -N8 -tu8 < /dev/urandom"));

▍ 4.4 Компиляция медленных функций


Существует утилита gp2c, с помощью которой можно скомпилировать функции в разделяемую библиотеку, которую можно потом подгрузить в gp.

По результатам моих экспериментов могу сказать, что это имеет смысл только в случаях, когда функция выполняемая параллельно большая, т. е. состоит из многих операторов. Тогда мы можем получить выигрыш в 3-10 раз. В среднем 3 раза. Но для маленьких функций выигрыш совсем небольшой. Около 25%. Поэтому к этому способу стоит обращаться в последнюю очередь.

И у него есть ряд ограничений. Вот они:

  1. Функция не должна использовать глобальные переменные.
  2. Экспорт переменных export(var1, ...) в скомпилированную функцию не работает.
  3. Не нужно использовать параллелизацию (parfor и т.п.) в компилируемой функции. Это сложно и нетривиально. У меня не вышло.
Итак, нужную для компиляции нам функцию (BruteForcePart) выносим в отдельный файл brute_force.gp и компилируем так:

gp2c -g -s _c brute_force.gp

После компиляции мы получаем .so файл, а также .c и .run файлы. Вы удивитесь, насколько похож полученный Си-код на исходный GP-код.

Эта команда сгенерирует библиотеку с поддержкой сборки мусора, и имя функции будет иметь суффикс _c, что даст нам возможность в одной и той же GP-программе сравнить скорость скомпилированной и обычной версии.

В скрипт нужно вставить содержимое файла .run. Оно подгружает в gp скомпилированную библиотеку. Выглядит примерно так:

install("init_brute_force","v","init_brute_force_c","./brute_force.gp.so");
install("BruteForcePart","D0,G,D0,G,D0,G,p","BruteForcePart_c","./brute_force.gp.so");

Хочу заметить, что компиляция в нативный код — это последнее средство. Так как могут всплыть такие тонкости, для которых потребуется изучать внутренности PARI. И порог входа сразу сильно взлетит. Но, впрочем, если вы будете следовать моим рекомендациям, то может и быстро получиться.

Я провозился целый день, чтобы разобраться, но на маленькой функции получил выигрыш всего в 25% по скорости. Так что занимайтесь этим, если у вас всё отлажено и вам нужно запустить многодневные вычисления, где действительно экономия по времени важна.

▍ 4.5 Другие параллельные конструкции


В этой части я обзорно покажу какие ещё есть конструкции в GP для параллельной работы. Детали смотрите в руководстве пользователя.
  • parapply(f, x) — параллельное вычисление функции f над элементами x
  • pareval(x) — параллельное вычисление элементов x, где x массив замыканий. Каждое замыкание должно быть без параметров, но внутри него можно использовать функции с параметрами
  • parforeach(V, x, expr1, {r}, {expr2 }) — Вычисляет параллельно выражение expr1 с аргументом x, где x проходит через все элементы V
  • parforprime(p = a, {b}, expr1, {r}, {expr2 }) — то же, что и parfor(), но p принимает значения только простых чисел
  • parforprimestep(p = a, {b}, q, expr1, {r}, {expr2 }) — то же, что и parfor(), но p принимает значения простых чисел из арифметической прогрессии a + kq
  • parsum(i = a, b, expr) — параллельно вычисляет сумму выражений expr, где i проходит от a до b в произвольном порядке
  • parvector(N, i, expr ) — Возвращает массив как и vector(N,i,expr), но инициализирует его параллельно
  • parselect(f, A, {flag = 0}) — Выбирает параллельно элементы A с помощью выбирающей функции f
  • parforvec(X = v, expr1, {j}, {expr2}, {flag}) — Вычисляет выражения expr2 (зависящие от X и j) для X создаваемого как в управляющей структуре forvec

5. Хитрости, которые я использую в Linux


Часто при целочисленном моделировании я вставляю в критическую секцию некую фильтрацию, а потом вывод на экран через print() или printf(), а вывод самого скрипта перенаправляю в текстовый файл.

▍ 5.1 Утилита sort


Потом для поиска наилучших результатов использую сортировку по одной из колонок, или не нескольким колонкам.

Примерно так:

gp -q transpose.gp >> trans.csv

Две правых угловых скобки означают, что я добавляю результаты в конец файла. Далее с помощью утилиты sort я ищу наилучшие результаты.

cat trans.csv | sort -t, -k4,4 -n > sorted.csv

В данном случае я использую как разделитель запятую (в случае пробела ключ -t не нужен), сортирую по четвёртой колонке -k4,4, данные в четвёртой колонке сортируются как числа, а не как строки -n.
sort — отличная утилита! Она сортирует многогигабайтные файлы с сумасшедшей скоростью.

▍ 5.2 Утилита uniq


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

Их нужно удалить после сортировки. Для этого служит полезная простая утилита uniq. Она обнаруживает только смежные строки, поэтому перед её использованием сортировка обязательна.

cat trans.csv | uniq > uniq.csv

В принципе у sort есть ключ -u, который тоже решает эту проблему.

6. Итоги


Параллельные вычисления теперь доступны каждому. С языком GP и стандартными утилитами типа sort использовать их не сложнее программирования на Basic (трудно придумать язык проще). И большое количество ядер даёт новое качество, позволяя экономить наше время жизни.

Если у вас нет, как у меня (конфигурация описана тут), под столом Epyc c 32-ядрами и 256GB RAM, то вы можете арендовать мощный сервер с 64 или даже 128 ядрами в RUVDS. И даже на несколько часов.

Разработчики opensource-программы PARI уже почти 40 лет стараются сделать вычисления удобными, точными и быстрыми для людей, у которых программирование не является основной специальностью. Я не встречал ещё более простой схемы распараллеливания вычислений, чем в PARI.

Пробуйте и, может быть, вы разделите мой энтузиазм!

Telegram-канал со скидками, розыгрышами призов и новостями IT ?

Комментарии (18)


  1. AndreyDmitriev
    08.10.2024 12:58
    #27392114
    +1

    Что-то у меня под Windows собрать до работоспособного состояния пример extgcd.c, что там идёт в комплекте не получилось. Студии отчаянно не нравится заголовочник divll_pre.h, а gcc (minGW64) в общем собирает, и оно даже запускается, но крешится где-то в районе gp_read_stream(stdin). Я верю, что под Линуксом оно всё ОК, но вот Windows что-то не срастается.

    Я пробовал и по-простому "gcc extgcd.c -llibpari -o extgcd" и так, как интернет рекомендует "gcc extgcd.c -Wall -Wextra -pedantic -llibpari -std=c11 -g -o extgcd", но нет.


    1. inetstar Автор
      08.10.2024 12:58
      #27392190

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

      Проще всего скомпилировать в библиотеку самую тяжёлую часть и загружать из pari. И даже это нужно делать только тогда, когда вся логика уже безукоризненно отлажена, а вычисления реально занимают часы или дни.


      1. AndreyDmitriev
        08.10.2024 12:58
        #27392288
        +1

        Я так понял вы решили начать по хардкору, написав на Си с использованием libpari. 

        Это я просто "с конца" начал. Так то я в LabVIEW программирую, и хотелось бы (допустим) внешнюю математику подключить библиотекой (суть DLL), иначе pari останется "вещью в себе". Я, конечно, могу (вероятно) через командную строку с ней общаться, но хотелось бы именно библиотекой. Там ещё надо будет посмотреть как данные пробрасывать, но это дело техники - если я успешно соберу консольное приложение с libpari, то соберу и библиотеку. Там ещё могут настичь подводные камушки с параллельностью, так как LabVIEW "параллельна" сама по себе, и в некоторых редких случаях внешняя многопоточная библиотека может доставить головную боль, но и это решаемо. А "пристреливаться" и отлаживаться, я конечно не в библиотеке буду. У меня в общем практической задачки в руках нет, да и дополнительная зависимость - не всегда хорошо, тут просто чисто любознательный интерес для расширения кругозора.


    1. inetstar Автор
      08.10.2024 12:58
      #27392984

      CC = cc
INCDIR = /usr/include/pari
LIBDIR = /usr/lib64
CFLAGS = -O -I$(INCDIR) -L$(LIBDIR)
$(CC) $(CFLAGS) -o trans2c trans2c.c -lpari -lm

      Я нашёл такую рекомендацию где-то в сети. Правда, для Linux.


      1. AndreyDmitriev
        08.10.2024 12:58
        #27393804
        +1

        А, не, там что-то с исходником не так. Я дома ещё раз проверил на совсем простом примере, всё ОК.

        Вот код в test.c, я хз сколько там надо в pari_init() заказывать, надо доки почитать, пусть будет с запасом:

        #include "pari/pari.h"

int main()
{
	pari_init(8000000,500000);
	pari_printf("Pari test - begin\n");

	GEN ver = pari_version(), Major = gel(ver,1), minor = gel(ver,2), p = gel(ver,3);
	pari_printf("Pari Version M = %Ps m = %Ps p = %Ps\n", Major, minor, p);
    printf("pari-%ld.%ld.%ld\n", itos(Major), itos(minor), itos(p));

	GEN M, D;
	M = mathilbert(5);
	D = det(M);
	pari_printf("M = %Ps\nD = %Ps\n", M, D);
	
	pari_close();
	printf("Pari test - end\n");	
	return 0;
}

        Компилируем просто:

        gcc test.c -l pari -o test

        Вывод из test.exe:

        Pari test - begin
Pari Version M = 2 m = 17 p = 0
M = [1, 1/2, 1/3, 1/4, 1/5; 
	1/2, 1/3, 1/4, 1/5, 1/6; 
	1/3, 1/4, 1/5, 1/6, 1/7; 
	1/4, 1/5, 1/6, 1/7, 1/8; 
	1/5, 1/6, 1/7, 1/8, 1/9]
D = 1/266716800000
pari-2.17.0
Pari test - end

        И да, детерминант матрицы Гильберта 5x5 считает верно.


  1. Zara6502
    08.10.2024 12:58
    #27392118
    +1

    Параллельные вычисления теперь доступны каждому. С языком GP

    Прочитал, но не понял это утверждение, оно написано так, словно параллельные вычисления были недоступны. В статье вы показали метод parfor который полностью эквивалентен Parallel.For/Parallel.Foreach в C# которые доступны всем давно. И так получается что писать ту самую параллельность вам придётся самому - а это очень специфичная тема и "познаниями бейсика" вы тут не обойдётесь.

    Пока мне не удалось увидеть в PARI ничего интересного по сравнению с C#, и если последний будет легко изучать тем кто видел C/C++/Java, то PARI - это свой специфичный синтаксис.

    В общем как-то статья не убедила, я, если честно, думал вы расскажите о чем-то, действительно уровня бейсика, где мне не придётся думать о памяти, общих переменных и т.п. и я просто отдам код, а оно само выполнится параллельно. Но это же не так.


    1. inetstar Автор
      08.10.2024 12:58
      #27392180
      +1

      Параллельные вычисления в языках общего назначения - это сложно и если нужно за полчаса написать программу, то вы не сделаете этого никогда не запускав C# до этого. В PARI достаточно, чтобы каждая процедура была автономной для параллельности.

      А c PARI сделаете. А учитывая, что PARI поддерживает неограниченную точность и огромное количество математических функций, то задача ещё упрощается.


      1. Zara6502
        08.10.2024 12:58
        #27394478

        Параллельные вычисления в языках общего назначения - это сложно и если нужно за полчаса написать программу, то вы не сделаете этого никогда не запускав C# до этого. В PARI достаточно, чтобы каждая процедура была автономной для параллельности.

        и в чем разница между parfor и Parallel.For/Parallel.Foreach в C#? Где конкретно в статье пишется что в PARI всё просто?

        А c PARI сделаете

        ну пока я вижу чтобы начать пользоваться PARI его нужно точно так же изучать, в тогда разница?

        А учитывая, что PARI поддерживает неограниченную точность и огромное количество математических функций, то задача ещё упрощается.

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


    1. inetstar Автор
      08.10.2024 12:58
      #27392206

      Кстати, о самой фразе, которую вы процитировали. Это разные предложения.

      Язык GP просто облегчает задачу, когда нужно что-то по-быстрому посчитать.

      Я, например, могу написать параллельный код на Go или Python, но на GP для меня будет проще, если мне не нужны какие-то тонкости.


      1. Zara6502
        08.10.2024 12:58
        #27394506

        Кстати, о самой фразе, которую вы процитировали. Это разные предложения.

        И?

        Язык GP просто облегчает задачу, когда нужно что-то по-быстрому посчитать.

        Каким образом?

        Я, например, могу написать параллельный код на Go или Python, но на GP для меня будет проще, если мне не нужны какие-то тонкости.

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


    1. inetstar Автор
      08.10.2024 12:58
      #27392274

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

      Но в критической секции можно результаты функций в какой-то глобальный массив записывать.


      1. Zara6502
        08.10.2024 12:58
        #27394536

        Фактически, я часто просто вывожу в критической секции результаты параллельных функций в поток вывода или файл, а потом анализирую. Очень просто

        бла-бла-бла, это просто слова, вы статью же для чего-то писали, даже в конце вроде бы сослались на то - ради чего "Пробуйте и, может быть, вы разделите мой энтузиазм!".

        У вас это энтузиазм, но статья не передаёт должной информации, чтобы этот энтузиазм появился у других. Вы точно так же как и для C# озвучили проблемы параллельных вычислений, точно такие же как и в C#, предложили (свели) всю параллельность к parfor - так что в PARI есть такое особенное? Я вам озвучил с какой надеждой я зашел читать статью, но я не увидел этого. Вы мне предложили запихать цикл в параллельный код, что я и сейчас могу делать в C#, причем очень много задач невозможно решить параллельным циклом не передавая информацию между итерациями. Это огромная проблема параллельных вычислений в целом и PARI ее не решает, вы же сами об этом пишете.

        Никаких заморочек с общими переменными

        И где это в статье? У вас вся параллельность, которая указана в заголовке, уместилась в один маленький раздел и в один пример.

        Вообще общие переменные использую только для чтения, как константы.

        Я точно так же это делаю в C#, но это сильное ограничение, например при обработке графов, когда у вас меняется вес в одной итерации и его нужно учитывать в другой - это огромная проблема параллельности вообще. Как её решает PARI?

        Но в критической секции можно результаты функций в какой-то глобальный массив записывать.

        Я могу это же делать в C#, собираю параллельную обработку в глобальный массив, а потом опять считаю, так как я не могу сразу работать с параллельными данными во время исполнения. Так в чем тут выигрыш PARI? Вы это как-то объяснили?

        Повторю еще раз - parfor это точно такая же конструкция как и Parallel.For/Parallel.Foreach в C#, почему я должен бросить знакомый C#, изучать PARI? Вы же наверное собирались побудить людей к этому своей статьёй?


        1. inetstar Автор
          08.10.2024 12:58
          #27395806

          Я добавил раздел 0, где написано кому может понадобиться PARI. А также добавил ещё раздел 4.5, где обзорно пробежался по другим функциям параллельности в PARI.

          Если сравнивать с C#, то порог вхождения несопоставим с PARI и значительно выше у C#.

          C# - это язык общего назначения, а PARI - вещь заточенная для быстрого решения математических/вычислительных проблем.


          1. Zara6502
            08.10.2024 12:58
            #27395838
            +1

            спасибо, то что нужно.

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


        1. AndreyDmitriev
          08.10.2024 12:58
          #27396074
          +1

          почему я должен бросить знакомый C#, изучать PARI?

          Вы знаете, тут больше всё-таки от задачи зависит, которую вы решаете. Ну вот, к примеру, надо вам вычислить детерминанты нескольких матриц матриц Гильберта очень большого размера или там факториал или ещё что, причём где-то параллельно. Соответственно (если не писать всё руками в С# с нуля) придётся воспользоваться подходящей математической библиотекой. Если Pari имеет всю необходимую для решения задачки математику, то можно воспользоваться ею. И тут есть несколько вариантов. Можно пользоваться Pari просто как алгебраическим калькулятором, оставаясь в рамках предложенного скриптового языка gp, тогда да, parfor и всё остальное — наше всё, вызывать это дело надо будет из командной строки, гнать вычисления в файл, тут всё просто. Либо, если надо тесно интегрировать в C# (или другую среду), то можно написать всё на Си, используя libpari.dll (кстати, конвертера gp2c я под Windows не вижу, но есть исходники), ну и подключать как динамическую библиотеку. Можно враппер написать, вот, к примеру для питона я вижу cypari/cypari2. В этом случае нам parfor вообще никаким боком не нужен (да мы его из Си/C#/Python и не можем использовать, разве что через gp_embedded_init()), тем не менее для параллелизации вычислений есть несколько возможностей (судя по примерам).

          Можно воспользоваться тем, что предоставляет Pari, через mt_queue*, тут запускается три потока для параллельной обработки N1 N2 и М:

          #include <pari/pari.h>

GEN
Cworker(GEN d, long kind) { return kind? det(d): Z_factor(d); }

int
main(void)
{
  long i, taskid, pending;
  GEN M,N1,N2, in,out, done;
  struct pari_mt pt;
  entree ep = {"_worker",0,(void*)Cworker,20,"GL",""};
  /* initialize PARI, postponing parallelism initialization */
  pari_init_opts(8000000,500000, INIT_JMPm|INIT_SIGm|INIT_DFTm|INIT_noIMTm);
  pari_add_function(&ep); /* add Cworker function to gp */
  pari_mt_init(); /* ... THEN initialize parallelism */
  /* Create inputs and room for output in main PARI stack */
  N1 = addis(int2n(256), 1); /* 2^256 + 1 */
  N2 = subis(int2n(193), 1); /* 2^193 - 1 */
  M = mathilbert(80);
  in  = mkvec3(mkvec2(N1,gen_0), mkvec2(N2,gen_0), mkvec2(M,gen_1));
  out = cgetg(4,t_VEC);
  /* Initialize parallel evaluation of Cworker */
  mt_queue_start(&pt, strtofunction("_worker"));
  for (i = 1; i <= 3 || pending; i++)
  { /* submit job (in) and get result (out) */
    mt_queue_submit(&pt, i, i<=3? gel(in,i): NULL);
    done = mt_queue_get(&pt, &taskid, &pending);
    if (done) gel(out,taskid) = done;
  }
  mt_queue_end(&pt); /* end parallelism */
  output(out); pari_close(); return 0;
}

          Либо можно воспользоваться старым добрым openmp, наброшенным поверх pari_thread_*(), тогда это вот так будет выглядеть:

          #include <pari/pari.h> /* Include PARI headers */

#include <omp.h>       /* Include OpenMP headers */

#define MAXTHREADS 3  /* Max number of parallel threads */

int
main(void)
{
  GEN M,N1,N2, F1,F2,D;
  struct pari_thread pth[MAXTHREADS];
  int numth = omp_get_max_threads(), i;
  /* Initialise the main PARI stack and global objects (gen_0, etc.) */
  pari_init(8000000,500000);
  if (numth > MAXTHREADS) {
    numth = MAXTHREADS;
    omp_set_num_threads(numth);
  }
  /* Compute in the main PARI stack */
  N1 = addis(int2n(256), 1); /* 2^256 + 1 */
  N2 = subis(int2n(193), 1); /* 2^193 - 1 */
  M = mathilbert(80);
  /*Allocate pari thread structures */
  for (i = 1; i < numth; i++) pari_thread_alloc(&pth[i],8000000,NULL);
#pragma omp parallel
  {
    int this_th = omp_get_thread_num();
    if (this_th) (void)pari_thread_start(&pth[this_th]);
#pragma omp sections
    {
#pragma omp section
      {
        F1 = factor(N1);
      }
#pragma omp section
      {
        F2 = factor(N2);
      }
#pragma omp section
      {
        D = det(M);
      }
    } /* omp sections */
    if (this_th) pari_thread_close();
  } /* omp parallel */
  pari_printf("F1=%Ps\nF2=%Ps\nlog(D)=%Ps\n", F1, F2, glog(D,3));
  for (i = 1; i < numth; i++) pari_thread_free(&pth[i]);
  return 0;
}

          Ну и для любителей хардкора можно на классических POSIX потоках (эту простыню я под спойлер засуну):

          Скрытый текст
          #include <pari/pari.h> /* Include PARI headers */

#include <pthread.h>   /* Include POSIX threads headers */

void *
mydet(void *arg)
{
  GEN F, M;
  /* Set up thread stack and get thread parameter */
  M = pari_thread_start((struct pari_thread*) arg);
  F = QM_det(M);
  /* Free memory used by the thread */
  pari_thread_close();
  return (void*)F;
}

void *
myfactor(void *arg)  /* same principle */
{
  GEN F, N;
  N = pari_thread_start((struct pari_thread*) arg);
  F = factor(N);
  pari_thread_close();
  return (void*)F;
}

int
main(void)
{
  long prec = DEFAULTPREC;
  GEN M1,M2, N1,N2, F1,F2, D1,D2;
  pthread_t th1, th2, th3, th4; /* POSIX-thread variables */
  struct pari_thread pth1, pth2, pth3, pth4; /* pari thread variables */

  /* Initialise the main PARI stack and global objects (gen_0, etc.) */
  pari_init(32000000,500000);
  /* Compute in the main PARI stack */
  N1 = addis(int2n(256), 1); /* 2^256 + 1 */
  N2 = subis(int2n(193), 1); /* 2^193 - 1 */
  M1 = mathilbert(149);
  M2 = mathilbert(150);
  /* Allocate pari thread structures */
  pari_thread_alloc(&pth1,8000000,N1);
  pari_thread_alloc(&pth2,8000000,N2);
  pari_thread_alloc(&pth3,32000000,M1);
  pari_thread_alloc(&pth4,32000000,M2);
  /* pthread_create() and pthread_join() are standard POSIX-thread
   * functions to start and get the result of threads. */
  pthread_create(&th1,NULL, &myfactor, (void*)&pth1);
  pthread_create(&th2,NULL, &myfactor, (void*)&pth2);
  pthread_create(&th3,NULL, &mydet,    (void*)&pth3);
  pthread_create(&th4,NULL, &mydet,    (void*)&pth4); /* Start 4 threads */
  pthread_join(th1,(void*)&F1);
  pthread_join(th2,(void*)&F2);
  pthread_join(th3,(void*)&D1);
  pthread_join(th4,(void*)&D2); /* Wait for termination, get the results */
  pari_printf("F1=%Ps\nF2=%Ps\nlog(D1)=%Ps\nlog(D2)=%Ps\n",
              F1,F2, glog(D1,prec),glog(D2,prec));
  pari_thread_free(&pth1);
  pari_thread_free(&pth2);
  pari_thread_free(&pth3);
  pari_thread_free(&pth4); /* clean up */
  return 0;
}

          Ну то есть при использовании Pari в многопоточной среде и интеграции со своей средой многопоточность среды остаётся как есть, она заменяет все эти "parfor", то есть если вы пользуетесь С#, то и будете использовать те же привычные Parallel.For/Parallel.Foreach, ну или я если в LabVIEW, то буду использовать свои конструкции. Я не думаю, что параллелизация parfor какая-то особенная, чудес не бывает, но всегда можно запустить бенчмарк на одинаковых данных и проверить. Как-то так.


          1. Zara6502
            08.10.2024 12:58
            #27396404

            Вы знаете, тут больше всё-таки от задачи зависит, которую вы решаете

            нет, вопрос был к оформлению статьи и к вопросам, которые будут задаваться неизбежно, из статьи совсем не следует, что PARI именно такой как описываете вы. Если бы следовало, то вопросы бы не возникли.

            простите, остальное читать не стал. Мы с автором эту тему закрыли, все всем довольны.


          1. inetstar Автор
            08.10.2024 12:58
            #27396442
            +1

            Спасибо большое за примеры параллельного кода!

            Когда я компилировал с помощью gp2c, то увидел, что libpari имеет свои средства для параллелизации. А точнее итераторы для параллельного исполнения parfor_t. И функции parfor_init, parfor_stop. Но у меня не вышло их использовать без доработки.

            Что нагенерил gp2c
            // То, что должно работать параллельно
GEN BruteForcePart(GEN pub, GEN cnt, GEN nMaxLevel)
{
  pari_sp ltop = avma;
  setrand(gpextern("od -vAn -N8 -tu8 < /dev/urandom"));
  /* Моя секретная функция */
  avma = ltop;
  return gen_0;
}

// Для parfor_init нужно замыкание. Оно создаётся в виде функции
GEN anon_0(GEN i, GEN pub, long prec)
{
  pari_sp ltop = avma;
  GEN p1 = gen_0;
  p1 = BruteForcePart(pub, nMaxCalcInThread, nMaxX, prec);
  p1 = gerepilecopy(ltop, p1);
  return p1;
}

// В этой функции мы запускаем parfor в GP-скрипте
void transpose(GEN pub, GEN a)         /* void */
{
  pari_sp ltop = avma;
  GEN nMaxCount = gen_0;
  
  nMaxCount = gmulgs(default0("nbthreads", NULL), 100);
  {
    pari_sp btop = avma;
    while (1)
    {
      {
        parfor_t iter;    /* parfor */
        parfor_init(&iter, gen_1, nMaxCount, strtoclosure("anon_0", 1, pub));
        {
          pari_sp btop = avma;
          GEN i = gen_0, r = gen_0, p1 = gen_0;
          while ((p1 = parfor_next(&iter)))
          {
            r = gcopy(gel(p1, 2));
            i = gcopy(gel(p1, 1));
            if (gequal0(r))
              continue;
            printf0("%s %s %d %d %d %d %d %s\n", mkvecn(8, gel(a, 1), gel(a, 2), gel(r, 1), lift(gel(gel(r, 2), 1)), lift(gel(gel(r, 2), 2)), lift(gel(pub, 1)), lift(gel(pub, 2)), gel(a, 5)));
            parfor_stop(&iter);
            goto label1;
            if (gc_needed(btop, 1))
              gerepileall(btop, 3, &i, &r, &p1);
          }
        }
      }
      avma = btop;
    }
    label1:;
  }
  avma = ltop;
  return;
}

            Этот код у меня не заработал. Подозреваю, что дело в том, что замыкание должно быть создано с помощью динамической компиляции GP-кода. Я не разобрался, как можно скомпилированную Си-функцию подсунуть в parfor_init.


            1. AndreyDmitriev
              08.10.2024 12:58
              #27396718
              +1

              О, спасибо огромное, а то я собрался уже было WSL расчехлять. Теперь понятнее.

              Я, кстати, попробовал навскидку gp_embedded_init()/gp_embedded() и эти функции работают, что даёт возможность напрямую вызывать скрипт вообще откуда угодно. Это не самый эфффективный способ (грубо говоря там интерпретатор работает), но удобно. Офигенная игрушка, на самом деле, спасибо за наводку.