Hi everyone! As you may already know (if you don't then I advice you to [learn](http://mirror.codeforces.com/blog/entry/22175) it), in 2D geometry it is convenient to use complex numbers to handling points operations and rotations. Now I would like to tell you about similar construction, which allows you to work efficiently with 3D geometry, in particular to use vectors and linear operations on them, calculate many popular operations (like dot or cross products) and maintain rotations in 3D space. ↵
↵
[cut]↵
<br>↵
↵
Let $\vec r_1 = (x_1, y_1, z_1), \vec r_2 = (x_2, y_2, z_2)$ — are vectors in 3D space. Following notation from analytical geometry will be used:↵
↵
1. Dot product: $(\vec r_1, \vec r_2) = x_1 x_2 + y_1 y_2 + z_1 z_2 = |\vec r_1| \cdot |\vec r_2| \cos \varphi$, where $\relax \varphi$ — is the angle between $\vec r_1$ and $\vec r_2$ in their common plane, $|\vec r_1|$ — length of vector $\vec r_1$.↵
↵
2. Cross product: $[\vec r_1, \vec r_2] = \begin{vmatrix} \vec i & \vec j & \vec k \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \end{vmatrix} = \begin{pmatrix}y_1 z_2 - z_1 y_2 \\ z_1 x_2 - x_1 z_2 \\ x_1 y_2 - y_1 x_2 \end{pmatrix}$, where $\vec i, \vec j, \vec k$ — unit vectors, directed along $\relax x$, $\relax y$ and $\relax z$ axes. In 3D space cross product is a vector having length $|\vec r_1| \cdot |\vec r_2| \sin \varphi$ and directed orthogonally to their common plane in such way that if you'll watch on vectors from the endpoint of cross product, shortest rotation from $\vec r_1$ to $\vec r_2$ will be counter-clockwise. ↵
↵
So, quaternion — is [hypercomplex number](https://en.wikipedia.org/wiki/Hypercomplex_number), which can be represented as $\relax q = s+xi+yj+zk$, where $\relax s, x, y, z$ — are real numbers and $\relax i, j, k$ — are imaginary units. We can define multiplication operation on quaternions with following equality $\relax {i^2 = j^2 = k^2 = i j k = -1}$. The whole quaternionic multiplication table can be derived from this identity:↵
↵
$$\begin{bmatrix}↵
\times & {\bf 1} & {\bf i} & {\bf j} & {\bf k} \\↵
{\bf 1} & 1 & i & j & k \\↵
{\bf i} & i & -1 & k & -j \\↵
{\bf j} & j & -k & -1 & i \\↵
{\bf k} & k & j & -i & -1↵
\end{bmatrix}$$↵
↵
Assuming that quaternion multiplication is distributive over addition ($\relax a (b+c) = ab+ac, (b+c)a=ba+ca$), we can multiply the quaternions, "opening parenthesis". Note that quaternionic multiplication will be associative ($\relax (ab)c=a(bc)$), but in common case not commutative (there are $\relax a, b$ such that $\relax ab \neq ba$).↵
↵
Quaternions can be also represented as $\relax s+ \vec v$, where $\relax \vec v = xi + yj + zk$ — is vector of 3D space with unit vectors $\relax i, j, k$. Components $\relax s$ and $\vec v$ are called scalar and vector parts of quaternions. Assume we have quaternions $a+\vec b$ and $c+\vec d$. Then their multiplication can be represented as $(a+\vec b)(c+\vec d) = ac + a \vec d + c \vec b + \vec b \vec d$. Consider multiplication of two quaternions with zero scalar parts.↵
↵
$\relax (x_{1} i + y_{1} j + z_{1} k)(x_{2} i + y_{2} j + z_{2} k) =$$\relax -(x_{1} x_{2} + y_{1} y_{2} + z_{1} z_{2}) + i(y_{1} z_{2} - z_{1} y_{2}) + j(z_{1} x_{2} - x_{1} z_{2}) + k(x_{1} y_{2} - y_{1} x_{2})$. ↵
↵
Note that it can be shortly written as $\vec {r_1} \vec {r_2} = [\vec{r_1}, \vec{r_2}] - (\vec{r_1}, \vec{r_2})$. Thus we have following formula for quaternionic multiplication:↵
↵
$$(a+\vec b)(c + \vec d) = ac - (\vec b, \vec d) + a\vec d + c \vec b + [\vec b, \vec d]$$↵
↵
Finally note that since $\relax ij = k$, quaternion $\relax a+bi+cj+dk$ can be represented in form $\relax (a+bi)+(c+di)j=A+Bj$, where $\relax A, B$ — are complex numbers. Also note that since $\relax ij=-ji$, we have identity $\relax (a+bi)j = j(a-bi)$ or, introducing notation $\hat A = a-bi$: $Aj=j \hat A$. Thus multiplication of quaternions $\relax A+Bj$ and $\relax C+Dj$ can be represented as $\relax AC+ADj+BjC+BjDj=(AC-B \hat D)+(AD+B\hat C)j$. Number $\hat A = a-bi$, introduced earlier, is called conjugated to complex number $\relax A$. So we have another one formula for quaternion multiplication:↵
↵
$$(A+Bj)(C+Dj) = (AC - B \hat D) + (AD + B \hat C)j$$↵
↵
Let's show that any non-zero quaternion is invertible, i.e. for each quaternion $\relax q$ there is inverse $\relax q^{-1}$ such that $\relax qq^{-1}=1$. In this case we can divide quaternions by multiplicating with inverse element. Like in complex numers we can introduce for quaternion $q=a+\vec b$ conjugated quaternion $\hat q = a - \vec b$. From that formula we can see that $q \hat q = a^2 + (\vec b, \vec b)$. Thus we can introduce quaternionic norm (which is generalisation of length): $||q|| = q \hat q$ and inverse element $q^{-1} = \dfrac{\hat q}{||q||}$. Note that $\widehat{ab}=\hat b \hat a$. Indeed from multiplication formula we get $(c - \vec d)(a - \vec b) = ac - (\vec b, \vec d) - a\vec d - c \vec b - [\vec b, \vec d]$. From here we can see $\relax ||ab||=||a||\cdot||b||$ and $\relax (ab)^{-1}=b^{-1} a^{-1}$.↵
↵
Now we are ready to learn about the most useful property of quaternions which is handling rotations in 3D space. From now on we will consider vectors to be quaternions with zero scalar part. Let's introduce conjugation operation of vector $a$ by quaternion $g$ the result of which will be vector <sup><a name = "ref10" href="#ref11">[1]</a></sup> $\relax \vec b=g \vec a g^{-1}$. This identity is equivalent to the following: $\vec b g = g\vec a$. Let $g = s + \vec v$, then we can get $s \vec b - (\vec b, \vec v) + [\vec b, \vec v] = s \vec a - (\vec v, \vec a) - [\vec a, \vec v]$. Considering scalar and vector parts separately we will get:↵
↵
$$\begin{cases}↵
(\vec b, \vec v) = (\vec a, \vec v) \\↵
s\vec b + [\vec b, \vec v] = s\vec a + [\vec v, \vec a]↵
\end{cases}$$↵
↵
Note that because $||ab|| = ||a|| \cdot ||b||$ we can get that norms of vectors $\vec a$ and $\vec b$ are equal. First equation of the system means that orthogonal projections on axis $\vec v$ are equal for $\vec a$ and $\vec b$. Set of vectors with same norm and projection on some direction form a circumference around that direction. Thus we see that $\vec b$ can be obtained from $\vec a$ by its rotation around $\vec v$ on some angle. Let's find it. From now on we will say that rotation is clockwise or counter-clockwise depending on what it's like if we watch it from the endpoint of vector around which rotation is happening.↵
↵
Let's assume that $\relax ||g||=1$. If it's not true, we can divide $\relax g$ by square root of its norm, it will not affect $\vec b$. Now we can assume that $g = \cos \varphi + \vec i \sin \varphi$, where $||\vec i|| = 1$. Also using first equation and using notation $\vec a' = \vec a - \vec i(\vec a, \vec i)$ and $\vec b' = \vec b - \vec i(\vec b, \vec i)$, we can transform second equation in following way: $\vec b' \cos \varphi + [\vec b', \vec i] \sin \varphi = \vec a' \cos \varphi + [\vec i, \vec a'] \sin \varphi$. (On the other hand we subtracted $\vec i (\vec a, \vec i) \cos \varphi = \vec i (\vec b, \vec i) \cos \varphi$ from both sides, on the other one since $\relax [a+b, c]=[a, c] + [b, c]$ we can see that $[\vec i, \vec a'] = [\vec i, \vec a + \vec i (\vec a, \vec i)] = [\vec i, \vec a]$)↵
↵
Note that $\vec a'$ and $\vec b'$ are lying in the plane orthogonal to $\vec v$, thus, this is the rotation plane. Now we should see that $[\vec b', \vec i]$ is vector $\vec b'$ rotated $90^\circ$ around vector $\vec i$ clockwise. Thus left part of second equation is $\vec b'$, rotated around $\vec i$ clockwise on angle $\relax \varphi$, and the right one is $\vec a'$, rotated around $\vec i$ counter-clockwise on $\relax \varphi$. Thus we see that $\vec b$ is $\vec a$ rotated around $\vec i$ counter-clockwise on $\relax 2 \varphi$.↵
<hr>↵
To sum up: quaternion with unit norm $g = \cos \left(\dfrac{\varphi}{2}\right) + \vec i \sin \left(\dfrac \varphi 2\right)$ has such property that for each $\vec a$ expression $g \vec a g^{-1}$ specifies vector $\vec a$ rotated around $\vec i$ on $\relax \varphi$ counter-clockwise.↵
↵
<span style="font-size:8pt"><a name = "ref11" href="#ref10">[1]</a> $(s+\vec v) \vec a (s - \vec v) = (s+\vec v)((\vec a, \vec v) + s\vec a - [\vec a, \vec v])$. Scalar part of this product equals to $s(\vec a, \vec v) - s(\vec a, \vec v) + (\vec v, [\vec a, \vec v]) = 0$.</span>↵
↵
And now the implementation. Ideone: [#BBOx9H](http://ideone.com/BBOx9H)↵
↵
```↵
// To begin, we introduce the notation for the class and its elements, which we will use.↵
typedef double ftype;↵
typedef complex<ftype> point;↵
typedef complex<point> quater;↵
#define qA real()↵
#define qB imag()↵
#define qs qA.real()↵
#define qx qA.imag()↵
#define qy qB.real()↵
#define qz qB.imag()↵
↵
const ftype pi = acos(-1.);↵
const ftype eps = 1e-9;↵
↵
// Multiplication via formulas in complex terms.↵
quater operator * (quater a, quater b)↵
{↵
return {{a.qA * b.qA - a.qB * conj(b.qB)},↵
{a.qA * b.qB + a.qB * conj(b.qA)}};↵
}↵
↵
// Conjugated element↵
quater conj(quater a)↵
{↵
return {conj(a.qA), - a.qB};↵
}↵
↵
// Quaternionic norm.↵
ftype norm(quater a)↵
{↵
return (a * conj(a)).qs;↵
}↵
↵
// Length.↵
ftype abs(quater a)↵
{↵
return sqrt(norm(a));↵
}↵
↵
// Dividing via inverse element.↵
quater operator / (quater a, quater b)↵
{↵
return a * conj(b) / point(norm(b));↵
}↵
↵
// Quaternion from vector coordinates.↵
quater vec(ftype x, ftype y, ftype z)↵
{↵
return {{0, x}, {y, z}};↵
}↵
↵
// Basis vectors.↵
const quater ex = vec(1, 0, 0);↵
const quater ey = vec(0, 1, 0);↵
const quater ez = vec(0, 0, 1);↵
↵
// Vector part of quaternion.↵
quater vec(quater a)↵
{↵
return a -= a.qs;↵
}↵
↵
// Dot product.↵
ftype dot(quater a, quater b)↵
{↵
return -(a * b).qs;↵
}↵
↵
// Cross product.↵
quater cross(quater a, quater b)↵
{↵
return vec(a * b);↵
}↵
↵
// Triple (mixed) product.↵
ftype mix(quater a, quater b, quater c)↵
{↵
return dot(a, cross(b, c));↵
}↵
↵
// Conjugation of vector a by quaternion g.↵
quater conj(quater a, quater g)↵
{↵
return g * a / g;↵
}↵
↵
// Quaternion representing rotation around i on phi counter-clockwise.↵
quater rotation(quater i, ftype phi)↵
{↵
return point(cos(phi / 2)) + i * point(sin(phi / 2));↵
}↵
↵
// Rotate a around i counter-clockwise.↵
quater rotate(quater a, quater i, ftype phi)↵
{↵
return conj(a, rotation(i, phi));↵
}↵
↵
// Check if two vectors are equal.↵
bool cmp(const quater &a, const quater &b)↵
{↵
return abs(a - b) < eps;↵
}↵
↵
// Any vector which is orthogonal to v.↵
quater ortho(quater v)↵
{↵
if(abs(v.qy) > eps)↵
return vec(v.qy, -v.qx, 0);↵
else↵
return vec(v.qz, 0, -v.qx);↵
}↵
↵
// Quaternion representing rotation from a to b via minimal angle.↵
quater min_rotation(quater a, quater b)↵
{↵
a /= abs(a);↵
b /= abs(b);↵
if(cmp(a, -b)) // Degenerate case :(↵
return rotation(ortho(b), pi);↵
else↵
return conj(a * (a + b));↵
}↵
↵
// Angle of rotation from [-pi; pi] defined by quaternion.↵
ftype get_angle(quater a)↵
{↵
a /= abs(a);↵
return remainder(2 * acos(a.qs), 2 * pi);↵
}↵
↵
// Quaternion representing rotation which turns nx in x, ny in y and [nx, ny] in z.↵
quater basis_rotation(quater nx, quater ny)↵
{↵
nx /= abs(nx);↵
ny /= abs(ny);↵
quater a = min_rotation(nx, ex);↵
ny = conj(ny, a);↵
quater b = min_rotation(ny, ey);↵
if(cmp(ny, -ey))↵
b = rotation(ex, pi);↵
return b * a;↵
}↵
```↵
↵
[cut]↵
<br>↵
↵
Let $\vec r_1 = (x_1, y_1, z_1), \vec r_2 = (x_2, y_2, z_2)$ — are vectors in 3D space. Following notation from analytical geometry will be used:↵
↵
1. Dot product: $(\vec r_1, \vec r_2) = x_1 x_2 + y_1 y_2 + z_1 z_2 = |\vec r_1| \cdot |\vec r_2| \cos \varphi$, where $\relax \varphi$ — is the angle between $\vec r_1$ and $\vec r_2$ in their common plane, $|\vec r_1|$ — length of vector $\vec r_1$.↵
↵
2. Cross product: $[\vec r_1, \vec r_2] = \begin{vmatrix} \vec i & \vec j & \vec k \\ x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \end{vmatrix} = \begin{pmatrix}y_1 z_2 - z_1 y_2 \\ z_1 x_2 - x_1 z_2 \\ x_1 y_2 - y_1 x_2 \end{pmatrix}$, where $\vec i, \vec j, \vec k$ — unit vectors, directed along $\relax x$, $\relax y$ and $\relax z$ axes. In 3D space cross product is a vector having length $|\vec r_1| \cdot |\vec r_2| \sin \varphi$ and directed orthogonally to their common plane in such way that if you'll watch on vectors from the endpoint of cross product, shortest rotation from $\vec r_1$ to $\vec r_2$ will be counter-clockwise. ↵
↵
So, quaternion — is [hypercomplex number](https://en.wikipedia.org/wiki/Hypercomplex_number), which can be represented as $\relax q = s+xi+yj+zk$, where $\relax s, x, y, z$ — are real numbers and $\relax i, j, k$ — are imaginary units. We can define multiplication operation on quaternions with following equality $\relax {i^2 = j^2 = k^2 = i j k = -1}$. The whole quaternionic multiplication table can be derived from this identity:↵
↵
\times & {\bf 1} & {\bf i} & {\bf j} & {\bf k} \\↵
{\bf 1} & 1 & i & j & k \\↵
{\bf i} & i & -1 & k & -j \\↵
{\bf j} & j & -k & -1 & i \\↵
{\bf k} & k & j & -i & -1↵
\end{bmatrix}$$↵
↵
Assuming that quaternion multiplication is distributive over addition ($\relax a (b+c) = ab+ac, (b+c)a=ba+ca$), we can multiply the quaternions, "opening parenthesis". Note that quaternionic multiplication will be associative ($\relax (ab)c=a(bc)$), but in common case not commutative (there are $\relax a, b$ such that $\relax ab \neq ba$).↵
↵
Quaternions can be also represented as $\relax s+ \vec v$, where $\relax \vec v = xi + yj + zk$ — is vector of 3D space with unit vectors $\relax i, j, k$. Components $\relax s$ and $\vec v$ are called scalar and vector parts of quaternions. Assume we have quaternions $a+\vec b$ and $c+\vec d$. Then their multiplication can be represented as $(a+\vec b)(c+\vec d) = ac + a \vec d + c \vec b + \vec b \vec d$. Consider multiplication of two quaternions with zero scalar parts.↵
↵
$\relax (x_{1} i + y_{1} j + z_{1} k)(x_{2} i + y_{2} j + z_{2} k) =$$\relax -(x_{1} x_{2} + y_{1} y_{2} + z_{1} z_{2}) + i(y_{1} z_{2} - z_{1} y_{2}) + j(z_{1} x_{2} - x_{1} z_{2}) + k(x_{1} y_{2} - y_{1} x_{2})$. ↵
↵
Note that it can be shortly written as $\vec {r_1} \vec {r_2} = [\vec{r_1}, \vec{r_2}] - (\vec{r_1}, \vec{r_2})$. Thus we have following formula for quaternionic multiplication:↵
↵
$$(a+\vec b)(c + \vec d) = ac - (\vec b, \vec d) + a\vec d + c \vec b + [\vec b, \vec d]$$↵
↵
Finally note that since $\relax ij = k$, quaternion $\relax a+bi+cj+dk$ can be represented in form $\relax (a+bi)+(c+di)j=A+Bj$, where $\relax A, B$ — are complex numbers. Also note that since $\relax ij=-ji$, we have identity $\relax (a+bi)j = j(a-bi)$ or, introducing notation $\hat A = a-bi$: $Aj=j \hat A$. Thus multiplication of quaternions $\relax A+Bj$ and $\relax C+Dj$ can be represented as $\relax AC+ADj+BjC+BjDj=(AC-B \hat D)+(AD+B\hat C)j$. Number $\hat A = a-bi$, introduced earlier, is called conjugated to complex number $\relax A$. So we have another one formula for quaternion multiplication:↵
↵
$$(A+Bj)(C+Dj) = (AC - B \hat D) + (AD + B \hat C)j$$↵
↵
Let's show that any non-zero quaternion is invertible, i.e. for each quaternion $\relax q$ there is inverse $\relax q^{-1}$ such that $\relax qq^{-1}=1$. In this case we can divide quaternions by multiplicating with inverse element. Like in complex numers we can introduce for quaternion $q=a+\vec b$ conjugated quaternion $\hat q = a - \vec b$. From that formula we can see that $q \hat q = a^2 + (\vec b, \vec b)$. Thus we can introduce quaternionic norm (which is generalisation of length): $||q|| = q \hat q$ and inverse element $q^{-1} = \dfrac{\hat q}{||q||}$. Note that $\widehat{ab}=\hat b \hat a$. Indeed from multiplication formula we get $(c - \vec d)(a - \vec b) = ac - (\vec b, \vec d) - a\vec d - c \vec b - [\vec b, \vec d]$. From here we can see $\relax ||ab||=||a||\cdot||b||$ and $\relax (ab)^{-1}=b^{-1} a^{-1}$.↵
↵
Now we are ready to learn about the most useful property of quaternions which is handling rotations in 3D space. From now on we will consider vectors to be quaternions with zero scalar part. Let's introduce conjugation operation of vector $a$ by quaternion $g$ the result of which will be vector <sup><a name = "ref10" href="#ref11">[1]</a></sup> $\relax \vec b=g \vec a g^{-1}$. This identity is equivalent to the following: $\vec b g = g\vec a$. Let $g = s + \vec v$, then we can get $s \vec b - (\vec b, \vec v) + [\vec b, \vec v] = s \vec a - (\vec v, \vec a) - [\vec a, \vec v]$. Considering scalar and vector parts separately we will get:↵
↵
$$\begin{cases}↵
(\vec b, \vec v) = (\vec a, \vec v) \\↵
s\vec b + [\vec b, \vec v] = s\vec a + [\vec v, \vec a]↵
\end{cases}$$↵
↵
Note that because $||ab|| = ||a|| \cdot ||b||$ we can get that norms of vectors $\vec a$ and $\vec b$ are equal. First equation of the system means that orthogonal projections on axis $\vec v$ are equal for $\vec a$ and $\vec b$. Set of vectors with same norm and projection on some direction form a circumference around that direction. Thus we see that $\vec b$ can be obtained from $\vec a$ by its rotation around $\vec v$ on some angle. Let's find it. From now on we will say that rotation is clockwise or counter-clockwise depending on what it's like if we watch it from the endpoint of vector around which rotation is happening.↵
↵
Let's assume that $\relax ||g||=1$. If it's not true, we can divide $\relax g$ by square root of its norm, it will not affect $\vec b$. Now we can assume that $g = \cos \varphi + \vec i \sin \varphi$, where $||\vec i|| = 1$. Also using first equation and using notation $\vec a' = \vec a - \vec i(\vec a, \vec i)$ and $\vec b' = \vec b - \vec i(\vec b, \vec i)$, we can transform second equation in following way: $\vec b' \cos \varphi + [\vec b', \vec i] \sin \varphi = \vec a' \cos \varphi + [\vec i, \vec a'] \sin \varphi$. (On the other hand we subtracted $\vec i (\vec a, \vec i) \cos \varphi = \vec i (\vec b, \vec i) \cos \varphi$ from both sides, on the other one since $\relax [a+b, c]=[a, c] + [b, c]$ we can see that $[\vec i, \vec a'] = [\vec i, \vec a + \vec i (\vec a, \vec i)] = [\vec i, \vec a]$)↵
↵
Note that $\vec a'$ and $\vec b'$ are lying in the plane orthogonal to $\vec v$, thus, this is the rotation plane. Now we should see that $[\vec b', \vec i]$ is vector $\vec b'$ rotated $90^\circ$ around vector $\vec i$ clockwise. Thus left part of second equation is $\vec b'$, rotated around $\vec i$ clockwise on angle $\relax \varphi$, and the right one is $\vec a'$, rotated around $\vec i$ counter-clockwise on $\relax \varphi$. Thus we see that $\vec b$ is $\vec a$ rotated around $\vec i$ counter-clockwise on $\relax 2 \varphi$.↵
<hr>↵
To sum up: quaternion with unit norm $g = \cos \left(\dfrac{\varphi}{2}\right) + \vec i \sin \left(\dfrac \varphi 2\right)$ has such property that for each $\vec a$ expression $g \vec a g^{-1}$ specifies vector $\vec a$ rotated around $\vec i$ on $\relax \varphi$ counter-clockwise.↵
↵
<span style="font-size:8pt"><a name = "ref11" href="#ref10">[1]</a> $(s+\vec v) \vec a (s - \vec v) = (s+\vec v)((\vec a, \vec v) + s\vec a - [\vec a, \vec v])$. Scalar part of this product equals to $s(\vec a, \vec v) - s(\vec a, \vec v) + (\vec v, [\vec a, \vec v]) = 0$.</span>↵
↵
And now the implementation. Ideone: [#BBOx9H](http://ideone.com/BBOx9H)↵
↵
```↵
// To begin, we introduce the notation for the class and its elements, which we will use.↵
typedef double ftype;↵
typedef complex<ftype> point;↵
typedef complex<point> quater;↵
#define qA real()↵
#define qB imag()↵
#define qs qA.real()↵
#define qx qA.imag()↵
#define qy qB.real()↵
#define qz qB.imag()↵
↵
const ftype pi = acos(-1.);↵
const ftype eps = 1e-9;↵
↵
// Multiplication via formulas in complex terms.↵
quater operator * (quater a, quater b)↵
{↵
return {{a.qA * b.qA - a.qB * conj(b.qB)},↵
{a.qA * b.qB + a.qB * conj(b.qA)}};↵
}↵
↵
// Conjugated element↵
quater conj(quater a)↵
{↵
return {conj(a.qA), - a.qB};↵
}↵
↵
// Quaternionic norm.↵
ftype norm(quater a)↵
{↵
return (a * conj(a)).qs;↵
}↵
↵
// Length.↵
ftype abs(quater a)↵
{↵
return sqrt(norm(a));↵
}↵
↵
// Dividing via inverse element.↵
quater operator / (quater a, quater b)↵
{↵
return a * conj(b) / point(norm(b));↵
}↵
↵
// Quaternion from vector coordinates.↵
quater vec(ftype x, ftype y, ftype z)↵
{↵
return {{0, x}, {y, z}};↵
}↵
↵
// Basis vectors.↵
const quater ex = vec(1, 0, 0);↵
const quater ey = vec(0, 1, 0);↵
const quater ez = vec(0, 0, 1);↵
↵
// Vector part of quaternion.↵
quater vec(quater a)↵
{↵
return a -= a.qs;↵
}↵
↵
// Dot product.↵
ftype dot(quater a, quater b)↵
{↵
return -(a * b).qs;↵
}↵
↵
// Cross product.↵
quater cross(quater a, quater b)↵
{↵
return vec(a * b);↵
}↵
↵
// Triple (mixed) product.↵
ftype mix(quater a, quater b, quater c)↵
{↵
return dot(a, cross(b, c));↵
}↵
↵
// Conjugation of vector a by quaternion g.↵
quater conj(quater a, quater g)↵
{↵
return g * a / g;↵
}↵
↵
// Quaternion representing rotation around i on phi counter-clockwise.↵
quater rotation(quater i, ftype phi)↵
{↵
return point(cos(phi / 2)) + i * point(sin(phi / 2));↵
}↵
↵
// Rotate a around i counter-clockwise.↵
quater rotate(quater a, quater i, ftype phi)↵
{↵
return conj(a, rotation(i, phi));↵
}↵
↵
// Check if two vectors are equal.↵
bool cmp(const quater &a, const quater &b)↵
{↵
return abs(a - b) < eps;↵
}↵
↵
// Any vector which is orthogonal to v.↵
quater ortho(quater v)↵
{↵
if(abs(v.qy) > eps)↵
return vec(v.qy, -v.qx, 0);↵
else↵
return vec(v.qz, 0, -v.qx);↵
}↵
↵
// Quaternion representing rotation from a to b via minimal angle.↵
quater min_rotation(quater a, quater b)↵
{↵
a /= abs(a);↵
b /= abs(b);↵
if(cmp(a, -b)) // Degenerate case :(↵
return rotation(ortho(b), pi);↵
else↵
return conj(a * (a + b));↵
}↵
↵
// Angle of rotation from [-pi; pi] defined by quaternion.↵
ftype get_angle(quater a)↵
{↵
a /= abs(a);↵
return remainder(2 * acos(a.qs), 2 * pi);↵
}↵
↵
// Quaternion representing rotation which turns nx in x, ny in y and [nx, ny] in z.↵
quater basis_rotation(quater nx, quater ny)↵
{↵
nx /= abs(nx);↵
ny /= abs(ny);↵
quater a = min_rotation(nx, ex);↵
ny = conj(ny, a);↵
quater b = min_rotation(ny, ey);↵
if(cmp(ny, -ey))↵
b = rotation(ex, pi);↵
return b * a;↵
}↵
```↵
