PARI-GP: как посчитать что-то просто, точно и параллельно

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

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

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

Оглавление

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 Компиляция медленных функций
5. Хитрости, которые я использую в Linux
  5.1 Утилита sort
  5.2 Утилита uniq
6. Итоги

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% по скорости. Так что занимайтесь этим, если у вас всё отлажено и вам нужно запустить многодневные вычисления, где действительно экономия по времени важна.

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.

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