Симуляции жидкости в реальном времени
57 #математика,#нравится,#рекомендую

Всем привет!

Наткнулся на видео о том, как работает симуляции жидкости:

И там же есть полезная ссылка на книгу о симуляции жидкости в реальном времени: Real-Time Fluid Dynamics for Games. Сишного кода в книге на сто строк всего. Безо всякой зубодробительной математики. Любопытно, что книга понятнее, чем видео. Схороню здесь, чтобы не потерять. Давно интересовался этой темой.

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

Возможно, в будущем, напишу реализацию жидкости/дыма на JavaScript и выложу на сайте с описанием и исходниками.

Рисуем сферу с помощью двух треугольников
269 #javascript,#three.js,#webgl,#интерактив,#математика,#своими руками

Сфера в трехмерной графике обычно состоит из сотни-другой треугольников, при этом половина из них не видна человеку, поскольку их отсекает face culling и/или zbuffer.

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

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

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

(далее…)

Практическое применение B-сплайна (B-Spline) в программировании
2K #python,#линейная интерполяция,#математика

Полтора дня разбирался в работе B-сплайна и как это применяется на практике. Хочу поделиться с вами результатами моих изысканий.

Немного теории

Допустим, у вас есть некий массив точек размером n. Что вам потребуется для работы B-сплайна:

  • определить количество интерполируемых точек (от этого зависит скорость работы, гладкость получившейся кривой и размер формулы для вычисления, количество точек равно k+1)
  • избавится от рекурсии формулы для вычисления

Общая формула для расчета коэффициентов:

    \begin{equation*}N_{i,k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}N_{i,k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}N_{i+1,k-1}(x)\end{equation*}

Формула рекурсивная и заканчивается на k равном единице:

    \begin{equation*} N_{i,1}(x)=\left\{\begin{array}{cr}1&\text{if }t_i\leq x\textless t_{i+1}\\0&\mathrm{otherwise}\quad\end{array}\right.\end{equation*}

Пояснения к формуле:

  • если делитель получается равным или близким по значению к нулю, значит вся дробь равна нулю
  • параметр i — это индекс текущей точки
  • x — значение от i до i + 1 с нужным вам шагом (например, i равно 7, и нам нужно создать 3 точки для интерполяции, значит берем x равным 7.0, 7.5, 8.0)
  • t — неубывающий массив индексов

Вычисление элементов массива t:

    \begin{equation*} t=\left\{\begin{array}{ll}0,&\mathrm{if}\ i<k\\i-k+1,&\mathrm{if}\ k\leq{i}\leq{n}\\n-k+2,&\mathrm{if}\ i>n\end{array}\right. \end{equation*}

Практика

Мне нужно было интерполировать по четырем точкам, поэтому я заранее на бумажке вычислил формулу без рекурсии:

    \begin{equation*} N_{i,3}(x)=\left\{\begin{array}{ll}0,&\mathrm{if}\ x\textless t_{i}\\\frac{x-t_i}{t_{i+2}-t_i}*\frac{x-ti}{t_{i+1}-ti},&\mathrm{if}\ t_i\leq x\textless t_{i+1}\\\frac{x-ti}{t_{i+2}-ti}*\frac{t_{i+2}-x}{t_{i+2}-t_{i+1}}+\frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}*\frac{x-ti}{t_{i+2}-t_{i+1}},&\mathrm{if}\ t_{i+1}\leq x\textless t_{i+2}\\\frac{t_{i+3}-x}{t_{i+3}-t_{i+1}}*\frac{t_{i+3}-x}{t_{i+3}-t_{i+2}},&\mathrm{if}\ t_{i+2}\leq x\textless t_{i+3}\\0,&\mathrm{if}\ x\geq t_{i+3}\end{array}\right. \end{equation*}

Про бумажку я, конечно, соврал, мне было удобнее сделать функцию на питоне вот такого вида:

def div(a, b):
    try:
        return a / float(b)
    except ZeroDivisionError:
        return 0
 
def N(i, k, x, t):
    if 1 == k:
        if t[i] <= x < t[i + 1]:
            return 1.0
        return 0.0
    a = div(x - t[i], t[i + k - 1] - t[i]) * N(i, k - 1, x, t)
    b = div(t[i + k] - x, x[i + k] - x[i + 1]) * N(i + 1, k - 1, x, t)
    return a + b

После этого добавил условие k==2 и все возможные возвращаемые значения в зависимости от x, затем условие k==3, значения которого зависят от k==2. Может вам будет удобнее все это сделать на бумажке.

Наконец можно приступать к интерполяции. Формула такая:

    \begin{equation*} \sum\limits_{i}^{i+k}=\vec{p_{i}}N_{i,k}(x) \end{equation*}

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

Небольшой пример

У нас есть некий массив точек, нам надо интерполировать с пятой по шестую точки с шагом 0.25, при этом $k$ равно 3 (значит для интерполяции нужно взять четыре точки).

    \begin{equation*} p_0=p_5*N_{0,3}(0.0)+p_6*N_{1,3}(0.0)+p_7*N_{2,3}(0.0)+p_8*N_{3,3}(0.0) \end{equation*}

    \begin{equation*} p_1=p_5*N_{0,3}(0.25)+p_6*N_{1,3}(0.25)+p_7*N_{2,3}(0.25)+p_8*N_{3,3}(0.25) \end{equation*}

    \begin{equation*} p_2=p_5*N_{0,3}(0.5)+p_6*N_{1,3}(0.5)+p_7*N_{2,3}(0.5)+p_8*N_{3,3}(0.5) \end{equation*}

    \begin{equation*} p_3=p_5*N_{0,3}(1.0)+p_6*N_{1,3}(1.0)+p_7*N_{2,3}(1.0)+p_8*N_{3,3}(1.0) \end{equation*}

Итого четыре новых точки вместо двух.

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

Блочное умножение матриц для программистов
9K #математика,#матрица

Ребята, сегодня я узнал об очень интересном способе перемножения матриц и спешу поделиться с вами. Как известно, матрицы умножаются друг на друга строка на столбец или столбец на строку. Формула поэлементного умножения матриц размерностью 2 на 2 выглядит вот так:

\begin{bmatrix}A_{11}&A_{12} \\A_{21}&A_{22}\end{bmatrix}\begin{bmatrix}B_{11}&B_{12}\\B_{21}&B_{22}\end{bmatrix}=\begin{bmatrix}C_{11}&C_{12} \\C_{21}&B_{22}\end{bmatrix},

где
\begin{bmatrix}C_{11}&C_{12} \\C_{21}&B_{22}\end{bmatrix}=\begin{bmatrix}A_{11}*B_{11}+A_{12}*B_{21}&A_{11}*B_{12}+A_{12}*B_{22} \\A_{21}*B_{11}+A_{22}*B_{21}&A_{21}*B_{12}+A_{22}*B_{22}\end{bmatrix}

Оказывается, с помощью этой же формулы можно рекурсивно разложить умножение больших матриц. Возьмем для примера умножение двух матриц 4×4:

\begin{bmatrix}a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}\end{bmatrix}\begin{bmatrix}b_{11}&b_{12}&b_{13}&b_{14}\\b_{21}&b_{22}&b_{23}&b_{24}\\b_{31}&b_{32}&b_{33}&b_{34}\\b_{41}&b_{42}&b_{43}&b_{44}\end{bmatrix}=\begin{bmatrix}c_{11}&c_{12}&c_{13}&c_{14}\\c_{21}&c_{22}&c_{23}&c_{24}\\c_{31}&c_{32}&c_{33}&c_{34}\\c_{41}&c_{42}&c_{43}&c_{44}\end{bmatrix}

Каждую матрицу можно разбить на блоки 2×2. По четыре блока на каждую:

\begin{array}{|cc|cc|}a_{11}&a_{12}&a_{13}&a_{14}\\a_{21}&a_{22}&a_{23}&a_{24}\\\hline a_{31}&a_{32}&a_{33}&a_{34}\\a_{41}&a_{42}&a_{43}&a_{44}\end{array}\ \ \ \begin{array}{|cc|cc|}b_{11}&b_{12}&b_{13}&b_{14}\\b_{21}&b_{22}&b_{23}&b_{24}\\\hline b_{31}&b_{32}&b_{33}&b_{34}\\b_{41}&b_{42}&b_{43}&b_{44}\end{array}\ \ \ \begin{array}{|cc|cc|}c_{11}&c_{12}&c_{13}&c_{14}\\c_{21}&c_{22}&c_{23}&c_{24}\\\hline c_{31}&c_{32}&c_{33}&c_{34}\\c_{41}&c_{42}&c_{43}&c_{44}\end{array}

Затем эти блоки можно умножить по формуле, приведенной в начале.

Пример для первого блока:

C_{11}=\begin{bmatrix}c_{11}&c_{12}\\c_{21}&c_{22}\end{bmatrix}=\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}\begin{bmatrix}b_{11}&b_{12}\\b_{21}&b_{22}\end{bmatrix}+\begin{bmatrix}a_{13}&a_{14}\\a_{23}&a_{24}\end{bmatrix}\begin{bmatrix}b_{31}&b_{32}\\b_{41}&b_{42}\end{bmatrix}

Таким же чудесным образом можно разбить матрицу 16×16, 256×256, и вообще любую квадратную матрицу с количеством строк/столбцов равным степени двойки.

Зачем это делать, спросит вдумчивый читатель?

Допустим у вас есть метод, который очень быстро умножает матрицы 4×4 (есть такие операции в NEON, SIMD и прочих MMX) и вам надо ускорить умножение огромной матрицы 32×32.

Воспользуйтесь кодом из лекции MIT (смотреть с 44-й минуты):

и вызывайте вместо

C[0] += A[0] * B[0];

свой метод, а вместо проверки (n == 1) поставьте (n == 4). Кстати, количество умножений в этом случае будет равно 512, а поэлементно — 32768. Почувствуйте разницу!

Оптимизированная версия кода из лекции MIT:

public void Rec_Mult(double[] C, int offC, double[] A, int offA, double[] B, int offB, int n, int rowsize) {
    if (n == 1) {
        C[offC] += A[offA] * B[offB];
    }
    else {
        final int d11 = 0;
        final int d12 = n / 2;
        final int d21 = (n / 2) * rowsize;
        final int d22 = (n / 2) * (rowsize + 1);
 
        final int C11 = offC + d11;
        final int A11 = offA + d11;
        final int B11 = offB + d11;
 
        final int C12 = offC + d12;
        final int A12 = offA + d12;
        final int B12 = offB + d12;
 
        final int C21 = offC + d21;
        final int A21 = offA + d21;
        final int B21 = offB + d21;
 
        final int C22 = offC + d22;
        final int A22 = offA + d22;
        final int B22 = offB + d22;
 
        // C11 += A11 * B11
        Rec_Mult(C, C11, A, A11, B, B11, n / 2, rowsize);
        // C11 += A12 * B21
        Rec_Mult(C, C11, A, A12, B, B21, n / 2, rowsize);
 
        // C12 += A11 * B12
        Rec_Mult(C, C12, A, A11, B, B12, n / 2, rowsize);
        // C12 += A12 * B22
        Rec_Mult(C, C12, A, A12, B, B22, n / 2, rowsize);
 
        // C21 += A21 * B11
        Rec_Mult(C, C21, A, A21, B, B11, n / 2, rowsize);
        // C21 += A22 * B21
        Rec_Mult(C, C21, A, A22, B, B21, n / 2, rowsize);
 
        // C22 += A21 * B12
        Rec_Mult(C, C22, A, A21, B, B12, n / 2, rowsize);
        // C22 += A22 * B22
        Rec_Mult(C, C22, A, A22, B, B22, n / 2, rowsize);
    }
}
Пример работы на C/C++ для (n == 2)

#include <stdio.h>
#include <string.h>
 
 
void Rec_Mult(int *C, const int *A, const int *B, int n, int rowsize)
{
    if (n == 2)
    {
        const int d11 = 0;
        const int d12 = 1;
        const int d21 = rowsize;
        const int d22 = rowsize + 1;
 
        C[d11] += A[d11] * B[d11] + A[d12] * B[d21];
        C[d12] += A[d11] * B[d12] + A[d12] * B[d22];
        C[d21] += A[d21] * B[d11] + A[d22] * B[d21];
        C[d22] += A[d21] * B[d12] + A[d22] * B[d22];
    }
    else
    {
        const int d11 = 0;
        const int d12 = n / 2;
        const int d21 = (n / 2) * rowsize;
        const int d22 = (n / 2) * (rowsize + 1);
 
        // C11 += A11 * B11
        Rec_Mult(C + d11, A + d11, B + d11, n / 2, rowsize);
        // C11 += A12 * B21
        Rec_Mult(C + d11, A + d12, B + d21, n / 2, rowsize);
 
        // C12 += A11 * B12
        Rec_Mult(C + d12, A + d11, B + d12, n / 2, rowsize);
        // C12 += A12 * B22
        Rec_Mult(C + d12, A + d12, B + d22, n / 2, rowsize);
 
        // C21 += A21 * B11
        Rec_Mult(C + d21, A + d21, B + d11, n / 2, rowsize);
        // C21 += A22 * B21
        Rec_Mult(C + d21, A + d22, B + d21, n / 2, rowsize);
 
        // C22 += A21 * B12
        Rec_Mult(C + d22, A + d21, B + d12, n / 2, rowsize);
        // C22 += A22 * B22
        Rec_Mult(C + d22, A + d22, B + d22, n / 2, rowsize);
    }
}
 
 
#define ROW_COUNT 8
 
 
void printMatrix(const char *name, const int *mat)
{
    printf("%s:\n", name);
 
    for (int i = 0; i < ROW_COUNT; ++i)
    {
        for (int j = 0; j < ROW_COUNT; ++j)
        {
            printf("%4d", mat[i * ROW_COUNT + j]);
        }
        printf("\n");
    }
    printf("\n");
}
 
 
int main()
{
    const int matA[ROW_COUNT * ROW_COUNT] =
    {
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 0, 0, 4, 5, 6,
        1, 2, 3, 1, 1, 4, 5, 6,
    };
 
    const int matB[ROW_COUNT * ROW_COUNT] =
    {
        2, 2, 2, 2, 2, 2, 2, 2,
        3, 3, 3, 3, 3, 3, 3, 3,
        0, 0, 0, 0, 0, 0, 0, 0,
        4, 0, 0, 0, 0, 0, 0, 2,
        4, 0, 0, 0, 0, 0, 0, 2,
        0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0,
        1, 1, 1, 1, 1, 1, 1, 1,
    };
 
    int matC[ROW_COUNT * ROW_COUNT];
    memset(matC, 0, sizeof(matC));
 
    Rec_Mult(matC, matA, matB, ROW_COUNT, ROW_COUNT);
 
    printMatrix("Matrix A", matA);
    printMatrix("Matrix B", matB);
    printMatrix("Multiply", matC);
 
    return 0;
}

Если немножко изменить код, то можно умножать неквадратные матрицы, но с шириной/высотой равной степени двойки (например, 16×4 на 4×16).

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

И напоследок — в процессе поиска информации по этому методу не нашел НИ ОДНОГО сайта на русском языке, который объяснял бы подобный метод для профанов. Одни только хвалебные оды самому себе в стиле: «реализовал блочное умножение матриц, алгоритм и реализацию ищите сами, у меня теперь все стало быстрее». Молодцы, так держать! Зачем делиться информацией, правда, ведь, да?

Простыми словами о представлении форматов float32 и double64 в памяти компьютера
6K #double64,#float32,#ieee754,#математика

Форматы float32 и double64 на пальцах

Вкратце, идея довольна простая: исходное число необходимо привести к нормализованному виду 1.NNN2 в двоичной системе счисления с помощью битового сдвига, затем записать дробную часть этого числа в мантиссу, а количество сдвигов в экспоненту. И завершающим штрихом сохранить знак исходного числа.

Формат float32 имеет такой вид: 1s8e23m, где s — количество бит под знак, e — экспонента, m — мантисса. Для формата double64 вид такой — 1s11e52m.

Первый бит кодирует знак числа. Если это ноль, число положительно, в противном случае число отрицательное.

Затем идут восемь бит экспоненты. Если совсем простыми словами, то экспонента — это количество сдвигов запятой в исходном числе, представленном в двоичном виде, для получения нормализованного числа вида 1.NNN2. Чтобы получить количество сдвигов из этих восьми бит надо отнять 12710. Отрицательная экспонента — сдвиг влево, положительная — сдвиг вправо.

После экспоненты следуют 23 бита мантиссы. Это дробная часть числа 1.NNN2.

А теперь практическая часть!

Перевод десятичной дроби во float32

Попробуем с полученными знаниями закодировать число -2.62510.

В двоичном виде это число имеет вид 10.1012. Для получения нормализованного числа необходимо запятую сдвинуть влево на один разряд. Получим число 1.01012 и экспоненту равной 1.

Для сохранения экспоненты к ней надо прибавить 12710. Получится 12610 или 011111102 в двоичном виде.

Берем нормализованное число 1.01012 и выделяем мантиссу 01012. Для сохранения этого числа, которое занимает 4 бита, надо добавить нули справа до 23 бит. Получится число 010100000000000000000002.

Знак числа отрицательный, значит первый бит равен единице.

Итог: s=12 e=011111102 m=010100000000000000000002, с чем я вас и поздравляю.

Обратный перевод: из float32 в число

Разберем пример из википедии. Есть число во float32 0xC000000016. В двоичной системе это будет 110000000000000000000000000000002.

Разобъем его на компоненты: s=12 e=100000002 m=000000000000000000000002.

Мантисса равна нулю, но, как уже было сказано выше, сохраняется только дробная часть мантиссы, а единица отбрасывается. Значит мантисса равна 1.000000000000000000000002.

Экспонента равна 100000002 или 12810, отнимаем 12710 и получается, что экспонента равна единице.

Возьмем мантиссу и сдвинем точку вправо на эту единицу, получится 10.00000000000000000000002, это 210 в десятичной системе счисления.

Знак числа равен единице, значит исходное число отрицательное.

Решение: -210, что и требовалось доказать.

P.S. Маленькие циферки справа от числа обозначают систему счисления, если кто не знает. Пример: два в десятичной системе — 210, один‑ноль‑один в двоичной системе — 1012. Кстати, красным цветом выделен знак числа, зелёным — экспонента, а мантисса, соответственно, синим.

Полезные ссылки

Блог Евгения Жирнова