math

Тригонометрия в бенчмарках

  • Дата публикации: 15 декабря, 2024
  • /Метки: benchmark, d, ldc, math
  • /Комментарии: 0

«Неважно, что ты любишь больше:
косинус ли, синус ли…»

Тригонометрия — основа многих приложений, от компьютерной графики до научных симуляций. Все мы привыкли вызывать sin и cos, не задумываясь, как они реализованы. А реализации могут быть разные! Работая над математической библиотекой для dlib2, я провел интересное исследование — какая тригонометрия лучше? Конечно, есть функции из std.math, и в большинстве случаев подойдут именно они. Но не все так просто — все зависит от того, что именно вы разрабатываете.

Если вы собираете обычное приложение, то кажется, что беспокоиться не о чем. Но если вам, по тем или иным причинам, нельзя обращаться к Phobos? Тогда есть два основных пути — sin и cos из стандартной библиотеки C, либо кастомная реализация, если код собирается под голое железо (например, при создании ядра ОС или программировании встраиваемой электроники). Но если вы используете LDC, то ничто не мешает использовать интринсики LLVM — они, оказывается, работают быстрее, чем std.math!

Я провел ряд тестов для всех вариантов тригонометрии:

  • Тест на точность — вычисление синуса и косинуса для 200 аргументов от -π до +π. Замерялась максимальная погрешность — расхождение результата с std.math.sin и std.math.cos;
  • Тест на производительность — время вычисления синуса и косинуса 1000000 раз.

Во всех кейсах я использовал LDC 1.39.0 под Windows 10. Получилось следующее:

  • std.math.sin, std.math.cos:
    • Время выполнения: 4 мс
  • LLVM интринсики llvm_sin, llvm_cos:
    • Время выполнения: 2 мс
    • Точность: абсолютная (макс. погрешность для sin: 0, для cos: 0)
  • Функции sin, cos из стандартной библиотеки C:
    • Время выполнения: 21 мс
    • Точность: абсолютная (макс. погрешность для sin: 0, для cos: 0)
  • Моя кастомная реализация на таблицах:
    • Время выполнения: 33 мс
    • Точность: порядка 10-7 (макс. погрешность для sin: 2.97038e-07, для cos: 1.78188e-07)

Также я пробовал версию с ассемблерными вставками, но она получилась почему-то медленнее кастомной — видимо, при использовании инлайнового ассемблера компилятор не задействует какие-то оптимизации (а еще есть мнение, что x87 fsin, fcos на современных процессорах медленные сами по себе). Смысла в таком варианте реализации особо нет, так что я его не стал рассматривать для включения в библиотеку.

В итоге в dlib2 войдут четыре реализации с таким приоритетом:

  • Если используется LDC, то синус и косинус — это интринсики (то есть, кодогенератор сам выбирает оптимальную реализацию под нужную архитектуру);
  • Если используются другие компиляторы (DMD, GDC):
    • Если код компилируется с поддержкой Phobos, то используются функции из std.math;
    • Если код собирается в режиме version(NoPhobos), но не version(FreeStanding) (то есть, под Windows или Unix-подобную ОС), то используются функции рантайма C;
    • Если же идет компиляция в bare metal, то используется кастомная реализация на таблицах.

Обновления

dlib 1.3.0

Вышла новая версия dlib. В библиотеке появился новый пакет dlib.math.random с реализацией генератора псевдослучайных чисел на основе C-функции rand. Проделан ряд улучшений в математическом пакете: добавлена поддержка компилятора GDC в модуль dlib.math.sse, появилась новая функция интерполяции bezierQuadratic.

Количество скачиваний dlib в реестре DUB достигло 1400 в месяц — рекордный показатель за все время существования проекта!

Подготовка к релизу Dagon 0.16

Новая версия Dagon планируется к выпуску совсем скоро — на днях внес ряд багфиксов и улучшений в физику Newton (в частности, исправлен прыжок контроллера персонажа на плоских поверхностях), а также добавил встроенную функцию создания скриншотов — Application.takeScreenshot.

Интерполяция на основе сигмоиды

Для одного из шейдеров на GLSL мне потребовалась «умная» интерполяция цветов с возможностью изменять резкость перехода от одного значения к другому — от полностью плавного (линейного) до дискретного. В итоге получилась вот такая функция, которую я вывел на основе рациональной сигмоиды — может быть, кому-то пригодится:

float sigmoid(float x, float k)
{
    return (x + x * k - k * 0.5 - 0.5) / 
        (abs(x * k * 4.0 - k * 2.0) - k + 1.0) + 0.5;
}

При k = 0 функция обращается в линейную, при k = 1 — разрывается в точке 0,5. Вы можете увидеть, как это работает, при помощи интерактивного графика на Desmos: https://www.desmos.com/calculator/s0cwcrtzvs.

Результат этой функции передается в привычный mix – то есть, вместо mix(c1, c2, t) пишем mix(c1, c2, sigmoid(t, k)). Получится, например, такое:

(градиенты гамма-скорректированы)

dlib 0.9.0 beta

Вышла бета-версия коллекции библиотек dlib 0.9.0. Из основных нововведений стоит отметить новый модуль dlib.math.tensor — реализацию тензоров с поддержкой как статического, так и динамического выделения памяти (еще один шаг, приближающий dlib к NumPy и Matlab). Также значительно улучшен пакет dlib.image: появился экспорт в BMP и TGA, двумерный foreach для изображений, диапазоны для окон и произвольных прямоугольных регионов. Новый пакет dlib.network, как планируется, будет содержать независимую от Phobos поддержку сети и веб-функциональность (пока в нем есть только парсер URL).

А еще Atrium был успешно портирован под FreeBSD.

dlib 0.5

Не так давно состоялось очередное крупное обновление коллекции библиотек dlib — вышла версия 0.5, наиболее значительным нововведением которой стала поддержка ручного управления памятью (РУП). Но — обо всем по порядку…

  • Новый модуль dlib.core.memory предоставляет средства для ручного выделения и высвобождения динамической памяти, независимые от сборщика мусора и основанные на malloc/free. Имеется поддержка структур, классов и массивов. При использовании классов рекомендуется использовать интерфейс ManuallyAllocatable и перегружать метод free, который ответственен за удаление объекта — в противном случае корректное удаление в некоторых случаях не гарантировано (например, при доступе через интерфейс или родительский класс).
  • Началась работа по переводу всей dlib на РУП. Так, загрузчики изрбражений (PNG, JPEG, TGA, BMP) в новой версии полностью независимы от сборщика мусора. Для этого активно используется паттерн абстрактной фабрики, ответственный за создание изображений  в памяти. Кстати, в загрузчике PNG значительно улучшена поддержка индексированных изображений, для них добавлена поддержка альфа-канала.
  • Кроме того, на РУП переведены некоторые контейнеры из dlib.container — BST, ассоциативный массив. Реализован полностью ручной динамический массив (dlib.container.array).
  • Еще одна новинка — ООП для структур (dlib.core.oop). Это экспериментальный модуль, реализующий для структур прототипный стиль ООП с поддержкой множественного наследования и параметрического полиморфизма. Полностью заменить классы он, конечно, не может, но окажется весьма полезен, если нужно создавать объекты с наследованием в стеке. В будущем планируется переписать некоторые внутренние механизмы dlib с использованием этой легковесной объектной системы.
  • В пакете dlib.math появилась поддержка дуальных кватернионов. Это частный случай алгербы Клиффорда, обобщение кватернионов на поле дуальных чисел. Их можно использовать, например, для описания движения тел в кинематике — один дуальный кватернион охватывает и перенос, и вращение. Кстати, реализация обычных кватернионов через инкапсуляцию теперь совместима с векторами.
  • Изменения коснулись и пакета вычислительной геометрии. Усеченная пирамида (dlib.geometry.frustum) теперь задается с нормалями ограничивающих плоскостей, указывающими наружу пирамиды. Подвергся изменению API проверки пересечения Frustum с AABB. Исправлены ошибки в реализации AABB и плоскости.

dlib 0.3

Состоялся релиз коллекции библиотек dlib 0.3. Нововведения этой версии:

  • Добавлены абстрактные потоки ввода/вывода (dlib.core.stream), независимые от Phobos, а также интерфейс файловой системы (dlib.filesystem) с готовыми реализациями для POSIX и Windows — этот интерфейс можно использовать, например, для построения виртуальных ФС.
  • Добавлена начальная поддержка HDRI в dlib.image (реализация формата изображений с плавающей запятой в dlib.image.hdri). Кроме того, обеспечена поддержка распараллеливания обработки изображений (dlib.image.parallel), добавлена поддержка чтения форматов TGA и BMP. Чтение/запись графических форматов теперь основаны на потоках, поэтому имеется возможность загружать изображения, например, напрямую из архивов.
  • Элементы матриц (dlib.math.matrix) теперь располагаются по столбцам, а не по строкам. Это серьезно нарушило обратную совместимость, но если вы не используете внутренние данные матриц и пользуетесь только внешним API, то это изменение не должно повлечь никаких проблем.

Более полный чейнджлог, а также исходники релиза вы можете найти на GitHub:
https://github.com/gecko0307/dlib/releases/tag/v0.3.0

Рефакторинг матриц в dlib

  • Дата публикации: 11 октября, 2013
  • /Метки: d, dlib, math, matrix
  • /Комментарии: 0

На днях состоялось грандиозное обновление пакета линейной алгебры dlib.math. Изменения коснулись, главным образом, реализации матриц. Если раньше матрицы 2×2, 3×3 и 4×4 имели каждая отдельную независимую реализацию, то теперь все они являются частными случаями обобщенной квадратной матрицы Matrix!(T,N) (где T – тип элементов, N – размерность). Она содержит все необходимые общие методы для матриц любого размера (нахождение определителя, нахождение обратной матрицы, нахождение матрицы миноров и алгебраических дополнений и т.д.), оптимизированные, где это возможно, для размерностей 2, 3 и 4. Таким образом, нынешние специализации Matrix2x2f, Matrix3x3f и Matrix4x4f практически идентичны их прежним аналогам.

Новая реализация создана с учетом обратной совместимости, но все-таки есть несколько критичных изменений:

1. Больше нет шаблонов Matrix2x2!(T), Matrix3x3!(T), Matrix4x4!(T). Используйте вместо них Matrix!(T,2), Matrix!(T,3) и Matrix!(T,4). При этом псевдонимы на специализации типа Matrix2x2f и Matrix4x4d сохранены;

2. Нет доступа к элементам матриц 4×4 через поля m*, t* и h*. Возможен только доступ через поля a*. Это справедливо для матриц любого размера:

a11 a12 a13 a14 .. a1N
a21 a22 a23 a24 .. a2N
a31 a32 a33 a34 .. a3N
a41 a42 a43 a44 .. a4N
 :   :   :   :  .
aN1 aN2 aN3 aN4  ' aNN

2. Все аффинные преобразования (функции rotationMatrix, translationMatrix и др.) и утилитарные функции для матриц вынесены в отдельный модуль dlib.math.affine. Там же находятся функции right, up, forward, translation, scaling, которые раньше были опрелены как методы в Matrix4x4!(T). Благодаря UFCS, их и теперь можно использовать как методы – однако все они теперь представляют собой свойства только для чтения. Пока они определены только для Matrix!(T,4), но в будущем функции базиса (right, up, forwartd) будут доступны и для Matrix!(T,3).

3. В целях обратной совместимости сохраняются модули dlib.math.matrix2x2, dlib.math.matrix3x3, dlib.math.matrix4x4, но они помечены как deprecated. Вместо них импортируйте dlib.math.matrix (и dlib.math.affine, если вам нужны аффинные преобразования)

2. Не рекомендуется использовать identityMatrix3x3!(T) и identityMatrix4x4!(T). Единичные матрицы создаются при помощи статического метода identity: например, Matrix3x3f.identity.

3. Не рекомендуется трансформировать векторы методом transform. Вместо этого лучше использовать умножение вектора на матрицу: Vector3f(1, 2, 3) * myMatrix.

4. Любые матрицы можно создавать при помощи функции-фабрики matrixf, которая автоматически определяет размерность на основе входных данных:

auto m1 = matrixf(
    8, 3, 2, 0,
    4, 0, 2, 0,
    1, 3, 3, 0,
    0, 0, 3, 1
);

Это выражение создаст матрицу типа Matrix!(float,4) и присвоит ее переменной m1.

Убедительная просьба всем пользователям dlib сообщить мне (в Issues в репозитории на GitHub, либо на почту – gecko0307@gmail.com), если будут обнаружены какие-то несостыковки и баги, связанные с данным рефакторингом матриц.

WolframAlpha

Совершенно случайно набрел на интереснейший сервис — http://www.wolframalpha.com. Это гибрид поисковика, базы знаний и вычислительной системы, понимающий запросы на литературном английском.

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

32*56-756

…и конвертер величин:

100 kilograms in pounds

Если вам нужна система линейной алгебры, а Matlab/Octave под рукой нет, то сервис поможет и здесь. Например, можно решать уравнения:

solve x^2+2x-1=0

Или рисовать графики функций:

plot f(x)=x^2

Поддерживаются операции с матрицами:

matrix inverse {{8, 3, 2, 0},{4, 0, 2, 0},{1, 3, 3, 0},{0, 0, 3, 1}} decimal

Также у WolframAlpha можно спрашивать какие-нибудь научные факты, например, возраст Вселенной:

universe age in years

Или статистические данные, например, крупнейшие города США:

biggest cities of usa

Сервис умеет показывать погоду в указанном городе или регионе:

weather in moscow

Я еще не успел ознакомиться со всеми возможностями этой замечательной системы — уверен, там есть еще много интересного.

Журнал «FPS» №25

Вышел 25 номер электронного PDF-журнала «FPS», посвященного разработке игр, программированию, компьютерной графике и звуку.

Читайте в этом номере:

> SIGGRAPH 2013: новости с выставки
> VIGAMUS — музей видеоигр
> Моделирование волос в Blender
> Физический движок своими руками, часть II
> Дуальные числа и автоматическое дифференцирование
> Фильтрация изображений в dlib
> История Git
> Linux — для геймеров?
> Полезные консольные команды в Linux

Номер доступен для онлайн-чтения и загрузки на сервисе Issuu.com, Документах Google и Dropbox.

Последние новости по проекту вы можете узнать в публичной странице журнала в социальной сети Google+: http://gplus.to/fpsmag. Добавляйте нас в круги, оставляйте свои комментарии и отписывайтесь в нашем сообществе.

Дуальные числа и касательная к кривой Безье

Дуальные числа, поддержка которых не так давно появилась в dlib (dlib.math.dual), обладают замечательным свойством: с их помощью можно реализовать автоматическое дифференцирование функций. Если производить вычисления не над вещественными, а над дуальными числами, то в вещественной части результата получается значение самой функции в заданной точке, а в дуальной – значение ее производной.
При этом если оформить функцию в виде шаблона, она без лишних телодвижений расширяется до множества дуальных чисел. Следующий пример показывает дифференцирование простейшей квадратичной функции:

import std.stdio;
import dlib.math.dual;

T parabola(T)(T x)
{
    return x*x;
}

void main()
{
    float x = 1.0f;

    Dualf eval = parabola(Dualf(x, 1.0f));
    float value = eval.re;
    float deriv = eval.du;
    
    writeln(deriv);
}

При запуске программа выдаст значение производной – 2 для точки 1. Правильность результата нетрудно проверить, зная формулу производной степенной функции: если f(x) = xn, то f ‘(x) = nxn-1. Следовательно, если f(x) = x2, то f ‘(x) = 2x.

Теперь начинается самое интересное. Попробуем вместо скалярных величин взять векторные и дифференцировать функцию кривой Безье (dlib.geometry.bezier) для двумерного случая:

import dlib.math.dual;
import dlib.math.vector;
import dlib.geometry.bezier;

alias DualVector2f = Vector!(Dualf, 2);

void main()
{
    float t = 0.5f;
    
    DualVector2f eval = bezierCurveFunc2D(
        DualVector2f(Dualf(0.0f), Dualf(0.0f)),
        DualVector2f(Dualf(1.0f), Dualf(1.0f)),
        DualVector2f(Dualf(2.0f), Dualf(1.0f)),
        DualVector2f(Dualf(3.0f), Dualf(0.0f)),
        Dualf(t, 1.0f));
}

Результирующий вектор eval будет содержать в вещественной части точку на кривой, а в дуальной – вектор касательной к кривой в этой точке, который нам остается только нормировать:

Vector2f point = Vector2f(eval.x.re, eval.y.re);
Vector2f tangent = Vector2f(eval.x.du, eval.y.du).normalized;

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