Алгебра строится над пространством.

В качестве исходного пространства используется проективная плоскость, обозначаемая как: $RP^2$. Базис проективной плоскости составляют:

• e₁ и e₂ — «настоящие» (евклидовы) вектора
• e₀ — проективный вектор

Базис проективного пространства: {e₁, e₂, e₀}

• Чтобы построить именно геометрическую алгебру, в пространстве вводится метрика
• Геометрическая алгебра, построенная конкретно над проективной плоскостью, и есть PGA2D
• Элементы геометрической алгебры (в т.ч. и PGA2D) называются мультивекторами

Вводимая метрика для базиса $RP^2$: {+1, +1, 0}

Такая метрика означает, что: e₁e₁ = 1, e₂e₂ = 1, e₀e₀ = 0

Введиние метрики порождает алгебру, базис которой (в случае PGA2D) состоит из 8 элементов.

В моём подходе все формулы выстроены для следующего базиса:

Базис алгебры: {1, e₁, e₂, e₀, e₂₀, e₀₁, e₁₂, e₀₁₂}
Мультивектор : {s, nx, ny, d , px , py , pw , ps  }

Запись {s, nx, ny, …} обозначает мультивектор [mv = s + nx*e₁ + ny*e₂ + ...]

В других учебных материалах могут брать немного другие базисы:

• может бьть другой порядок базисных векторов, например: {1, e₀, e₁, e₂, ...}
• может быть переставлен порядок векторов в бивекторе/тривекторе: например, за базисный может браться не e₂₀, а e₀₂

Полная таблица перемножения базисных n-векторов:

Определения, использующиеся во всех видах GA:

• k-Vector — MV, состоящий только из Grade-k элементов
• k-Blade u (simple k-Vector u) — u можно разложить во внешнее произведение k штук 1-векторов: u = v₁ Λ v₂ Λ …
• Versor u — u можно разложить в геометрическое произведение k штук 1-векторов: u = v₁v₂…
• Все 1-векторы являются и блейдами и версорами
• Все блейды являются версорами

Применительно к PGA2D:

• Все E.Lines, I.Lines, E.Points, I.Points являются блейдами и версорами
• Motors в общем случае являются версорами, но не обязательно блейдами

Обоснование:

• E.Lines и I.Lines являются 1-векторами
• E.Points являются внешними произведениями линий: E.Point = E.Line1 Λ E.Line2
• I.Points:
  • можно разложить на произведение [px*e₂₀ + py*e₀₁ = e₀(-px*e₂ + py*e₁)]
  • исходя из свойств внешней алгебры (для любой GA), для базисных векторов геометрическое произведение заменяется внешним: e₀e₂ = e₀Λe₂ (аналогично для e₀e₁)
  • итого, мы успешно разложили I.Point на внешнее произведение 1-векторов
• Motors являются композициями отражений (т.е. геометрическими произведениями 1-векторов)

Все формулы ниже даются для базисов:

• Базис пространства: {e₁, e₂, e₀}
• Метрика:            {+1, +1, 0 }
• Базис алгебры:      {1, e₁, e₂, e₀, e₂₀, e₀₁, e₁₂, e₀₁₂}
• Мультивектор :      {s, nx, ny, d , px , py , pw , ps  }

Объекты алгебры PGA2D:

(E = Euclidean, I = Ideal)

Полезно ввести константы:

(US обозначает простой скаляр, но подчёркивает, что этот скаляр — часть алгебры)

Мультивекторы в формулах записываются в виде списков из 8 элементов:

 mv  = {s , nx , ny , d , px , py , pw , ps } ;
 mv1 = {s1, nx1, ny1, d1, px1, py1, pw1, ps1} ;
 mv2 = {s2, nx2, ny2, d2, px2, py2, pw2, ps2} ;

В формулах для произведений предполагаются:

• входные вектора: mv1, mv2
• результирующий вектор: mv

При перемножении нескольких векторов, например a Λ b Λ c :

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

Формула выводится напрямую из метрики:

• По определению GP:
  • eiej = —ejei
  • eiei = 0
• По метрике: e1e1 = +1, e2e2 = +1, e0e0 = 0

Далее просто подставляем и группируем:
(s1 + nx1∗e1 + ..)(s2 + nx2∗e2 + ..) =
(s1∗s2) + nx1∗s2∗(e1) + nx2∗s1∗(e2) + nx1∗nx2∗(e1e2) + ...

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

GP[mv1, mv2] := {s, nx, ny, d, px, py, pw, ps}

s  = nx1*nx2+ny1*ny2-pw1*pw2+s1*s2
nx = ny2*pw1-ny1*pw2+nx2*s1+nx1*s2
ny = -nx2*pw1+nx1*pw2+ny2*s1+ny1*s2
d  = -ps2*pw1-ps1*pw2-ny2*px1+ny1*px2+nx2*py1-nx1*py2+d2*s1+d1*s2
px = d2*ny1-d1*ny2+nx2*ps1+nx1*ps2-pw2*py1+pw1*py2+px2*s1+px1*s2
py = -d2*nx1+d1*nx2+ny2*ps1+ny1*ps2+pw2*px1-pw1*px2+py2*s1+py1*s2
pw = -nx2*ny1+nx1*ny2+pw2*s1+pw1*s2
ps = d2*pw1+d1*pw2+nx2*px1+nx1*px2+ny2*py1+ny1*py2+ps2*s1+ps1*s2

GP ассоциативно, т.е. (ab)c = a(bc)

Используется определение внешнего произведения на n-векторах a и b:

• $a \land b = \langle ab \rangle_{s+t}$
• здесь: Grade[a] = s, Grade[b] = t

Эта формула далее просто напрямую расширяется на произвольный мультивектор. И проводится общее суммирование по грейдам.

При этом если s+t ≥ 4, то произведение считается равным нулю.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

OP[mv1, mv2] := {s, nx, ny, d, px, py, pw, ps}

s  = s1*s2
nx = nx2*s1+nx1*s2
ny = ny2*s1+ny1*s2
d  = d2*s1+d1*s2
px = d2*ny1-d1*ny2+px2*s1+px1*s2
py = -d2*nx1+d1*nx2+py2*s1+py1*s2
pw = -nx2*ny1+nx1*ny2+pw2*s1+pw1*s2
ps = d2*pw1+d1*pw2+nx2*px1+nx1*px2+ny2*py1+ny1*py2+ps2*s1+ps1*s2

OP ассоциативно, т.е. (a Λ b) Λ c = a Λ (b Λ c)

Используется определение внутреннего произведения на n-векторах:

• $a \land b = \langle ab \rangle_{|s-t|}$
• здесь: Grade[a] = s, Grade[b] = t

Аналогично OP, эта формула далее расширяется на произвольный мультивектор.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

IP[mv1, mv2] := {s, nx, ny, d, px, py, pw, ps}

s  =  nx1*nx2+ny1*ny2-pw1*pw2+s1*s2
nx =  ny2*pw1-ny1*pw2+nx2*s1+nx1*s2
ny = -nx2*pw1+nx1*pw2+ny2*s1+ny1*s2
d  = -ps2*pw1-ps1*pw2-ny2*px1+ny1*px2+nx2*py1-nx1*py2+d2*s1+d1*s2
px =  nx2*ps1+nx1*ps2+px2*s1+px1*s2
py =  ny2*ps1+ny1*ps2+py2*s1+py1*s2
pw =  pw2*s1+pw1*s2
ps =  ps2*s1+ps1*s2

IP не ассоциативно, т.е. (a·b)·c ≠ a·(b·c)

Для произвольных мультивекторов a и b:

• $a \lor b = (a^{\ast} \land b^{\ast})^{\ast}$
• «*» здесь обозначает дуализацию (см. Dual)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

RP[mv1, mv2] := {s, nx, ny, d, px, py, pw, ps}

s  = d2*pw1+d1*pw2+nx2*px1+nx1*px2+ny2*py1+ny1*py2+ps2*s1+ps1*s2
nx = nx2*ps1+nx1*ps2+pw2*py1-pw1*py2
ny = ny2*ps1+ny1*ps2-pw2*px1+pw1*px2
d  = d2*ps1+d1*ps2-px2*py1+px1*py2
px = ps2*px1+ps1*px2
py = ps2*py1+ps1*py2
pw = ps2*pw1+ps1*pw2
ps = ps1*ps2

RP ассоциативно, т.е. (a V b) V c = a V (b V c)

Определение сэндвич-произведения для произвольных МВ:
WP[a, b] = bab~

(здесь «a~» это Reverse a, см. ниже)

Как следствие, для композиции моторов:
WP[a, m₂m₁] = m₂m₁am₁~m₂~

Вычислительно, рекомендую реализовывать функцию WP в виде: WP[a, {b1, b2, …}], где bi — список произвольной длины

Операция Grade projection просто берёт от МВ соответствующий грейд, зануляя остальные элементы.

Обозначается угловыми скобками. Например, обозначение Grade-2 проекции: $ \langle a \rangle_{2}$

1
2
3
4

Grade[mv, 0] := {s, 0 , 0 , 0, 0 , 0 , 0 , 0 } ;
Grade[mv, 1] := {0, nx, ny, d, 0 , 0 , 0 , 0 } ;
Grade[mv, 2] := {0, 0 , 0 , 0, px, py, pw, 0 } ;
Grade[mv, 3] := {0, 0 , 0 , 0, 0 , 0 , 0 , ps} ;

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

1
2

Dual[{s, nx, ny, d, px, py, pw, ps}] :=
    {ps, px, py, pw, nx, ny, d, s} ;

Реверс по определению «переворачивает» порядок произведения базисных векторов.

Например: (e₀e₁e₂)~ = e₂e₁e₀ = -e₂e₀e₁ = e₀e₂e₁ = -e₀e₁e₂

И аналогично для всех базисных элементов.

Конкретно для PGA2D, реверс просто меняет знак 2- и 3-грейд составляющих мультивектора.

1
2

Reverse[{s, nx, ny, d, px, py, pw, ps}] :=
    {s, nx, ny, d, -px, -py, -pw, -ps} ;

По определению, для произвольного мультивектора a:

• $a^{\perp} = GP[a, UPS]$
• Словами: полярный мультивектор суть геометрическое произведение исходного мультивектора на Unit Pseudoscalar
• UPS = {0, 0, 0, 0, 0, 0, 0, 1}

1
2

Polar[{s, nx, ny, d, px, py, pw, ps}] := 
    {0, 0, 0, -pw, nx, ny, 0, s}

В частности, поляризация E.Line порождает I.Point (Direction):

1
2

Polar[{0, nx, ny, d, 0 , 0 , 0, 0}] =
	{0, 0 , 0 , 0, nx, ny, 0, 0}

Про версоры:

• Версор — МВ, представимый в виде геом.произведения 1-векторов
• В PGA2D версорами являются: E.Line, I.Line, E.Point, I.Point, Motors

Вычисление нормы опирается на вычисление GP[mv~,mv]:

1
2
3
4
5
6

GP[Reverse[mv], mv] := {s, nx, ny, d, 0, 0, 0, 0}

s  = s*s + nx*nx + ny*ny + pw*pw
nx = -2*ny*pw + 2*nx*s
ny =  2*nx*pw + 2*ny*s
d  =  2*ps*pw + 2*ny*px - 2*nx*py + 2*d*s

Для указанных выше версоров GP[mv~,mv] имеет только скалярную составляющую, поэтому для версоров используют формулу $\langle GP[\tilde{vrsr}, vrsr] \rangle _{0}$ :

1

Grade[GP[Reverse[vrsr], vrsr], 0] := s*s + nx*nx + ny*ny + pw*pw

Обратный мультивектор по определению: GP[a, a⁻¹] = 1

Известна формула для версоров:
a⁻¹ = a~/GP[a~,a]

1

Inverse[mv] := Reverse[mv]/Norm[mv]^2

(здесь заложено, что Norm тоже работает только на версорах)

Строгое определение для произвольного МВ a:
$Norm[a] = \sqrt{GP[a, \tilde{a}]}$

Для версоров GP[a,a~] имеет только скалярную составляющую, т.е. для них можно использовать формулу:
$Norm[a] = \sqrt{\langle GP[a, \tilde{a}] \rangle _{0}}$

1

Norm[mv] := Sqrt[s*s + nx*nx + ny*ny + pw*pw]

Следствия:

• Norm[E.Line ] = Sqrt[nx² + ny²]
• Norm[E.Point] = |pw|
• Norm[I.Line ] = 0
• Norm[I.Point] = 0
• Norm[Rotator] = Sqrt[s² + pw²]
• Norm[Translator] = |s|

Строгое определение для произвольного МВ a:
$INorm[a] = \sqrt{GP[a^{\ast}, \tilde{a^{\ast}}]}$

Для версоров GP[a*,(a*)~] имеет только скалярную составляющую, т.е. для них:

1

INorm[mv] := Sqrt[d*d + px*px + py*py + ps*ps]

По определению, нормированный мультивектор удовлетворяет условию:
GP[a,a~] = 1

Поэтому, для версоров с Norm[a] ≠ 0 (т.е. для E.Point, E.Line и Motors):

1

Normalize[a] := a/Norm[a]

Для I.Point и I.Line вводится идеальная нормализация:

1

INormalize[a] := a/INorm[a]

Моторы вычисляются как экспоненты 2-векторов, для чего используется формула:

Exp[mv] = 1 + mv/1! + mv²/2! + mv³/3! + ...
(здесь возведение в степень — это геометрическое произведение вектора на себя)

А уж почему эспоненцирование 2-вектора является мотором — вопрос отдельный...

Формулы применения моторов (juxta означает GP):

• Композиция моторов m₁ и m₂:
  • Для нормированных/ненормированных моторов: m₂m₁
  • Если исходные моторы нормированные, то и результирующий будет нормированным
• Применить мотор m к объекту A: mAm~
  • Ненормированный мотор может менять норму/ид-норму A (проверить)
  • Ненормированный объект может менять норму/ид-норму при действии нормированных моторов (проверить)

В PGA2D существуют только 2 вида моторов (если не касаться флекторов): ротаторы и трансляторы. В общем случае композицией ротатора и транслятора является ротатор вокруг новой точки.

Смена знака всех элементов мотора не влияет на его действие

Ротатор на угол α вокруг нормированной E.Point NEP:

• NEP = {0, 0, 0, 0, px, py, 1, 0}, при этом д/б √(px² + py²) ≠ 0
• Rotator = Exp[α*NEP/2] = Cos[α/2] + Sin[α/2]*NEP
• Получившийся ротатор будет нормированным

Для ненормированной E.Point EP:

1
2
3

Rotator[α, EP] :=
    US*Cos[α/2] + Normalize[EP]*Sin[α/2]
(* получившийся ротатор будет нормированным *)

Транслятор на расстояние dist перпендикулярно ид-нормированной ид-точке I.Point INIP:

• INIP = {0, 0, 0, 0, px, py, 0, 0}, при этом √(px² + py²) = 1
• Exp[dist*INIP/2] = 1 + dist*INIP/2
• Получившийся транслятор будет нормированным

Для не-ид-нормированной ид-точки I.Point IP:

1
2

Translator[dist, IP] := US + (dist/2)*INormalize[IP]
(* получившийся транслятор будет нормированным *)

Кстати, I.Point можно получить через поляризацию E.Line или I.Line:

• Polar[{0, nx, ny, d, 0, 0, 0, 0}] = {0, 0, 0, 0, nx, ny, 0, 0}
• Нормированная E.Line породит ид-нормированную I.Point
• Ид-нормированная I.Line породит ид-нормированную I.Point

Деконструирование восстанавливает dist и ид-нормированную I.Point:
{0, 0, 0, 0, px, py, 0, 0}

Вывод формулы сделаем для общего случая ненормированного транслятора:

• Normed     translator: {1, 0, 0, 0, dist/2\*px, dist/2\*py, 0, 0}
• Non-normed translator: {A, 0, 0, 0, X, Y, 0, 0}

Деконструирование:

• В общем случае tr := {A, 0, 0, 0, X, Y, 0, 0}
• Если A<0, то просто работаем с tr := -1*tr
• Надо разрешить 3 уравнения:
  • X/A = px*dist/2
  • Y/A = py*dist/2
  • √(px² + py²) = 1

В итоге получаем:

• dist = 2*√((X/A)² + (Y/A)²)
• px = 2*(X/A)/dist
• py = 2*(Y/A)/dist

Деконструирование восстанавливает угол поворота alpha и нормированную E.Point:
{0, 0, 0, 0, px, py, 1, 0}

Вывод формулы сделаем для общего случая ненормированного ротатора:

• Normed     rotator: {Cos[α/2], 0, 0, 0, px*Sin[α/2], py*Sin[α/2], Sin[α/2], 0}
• Non-normed rotator: {C, 0, 0, 0, X, Y, S, 0}

Деконструирование:

• В общем случае rot := {C, 0, 0, 0, X, Y, S, 0}
• Надо разрешить 3 уравнения:
  • cc = Cos[α/2] = C/√(C² + S²)
  • ss = Sin[α/2] = S/√(C² + S²)
  • px*Sin[α/2] = X/√(C² + S²)
  • py*Sin[α/2] = Y/√(C² + S²)

В итоге получаем:

• Первые 2 уравнения содержат достаточно информации для извлечения α ∈ (-π,π] (например, в Python можно использовать α = 2*math.atan2(ss, cc))
• px = X/S
• py = Y/S

Ротатор {Cos[α/2], 0, 0, 0, px*Sin[α/2], py*Sin[α/2], Sin[α/2], 0}:

• При alpha = 0 ротатор переходит в Scalar
• При alpha = ±π ротатор переходит в E.Point

Транслятор {1, 0, 0, 0, dist/2*px, dist/2*py, 0, 0}:

• При dist = 0 транслятор переходит в Scalar

Пусть A и B — либо 2 E.Point, либо 2 E.Line.

Чтобы перевести A в B надо приложить мотор motor:

• B = WP[A, motor]
• motor = √(B/A) = Sqrt[GP[B, Inverse[A]]