Hi everyone! As you may already know (if you don't then I advice you to learn 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.

Let — are vectors in 3D space. Following notation from analytical geometry will be used:

1. Dot product: , where — is the angle between and in their common plane, — length of vector .

2. Cross product: , where — unit vectors, directed along , and axes. In 3D space cross product is a vector having length 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 to will be counter-clockwise.

So, quaternion — is hypercomplex number, which can be represented as , where — are real numbers and — are imaginary units. We can define multiplication operation on quaternions with following equality . The whole quaternionic multiplication table can be derived from this identity:

Assuming that quaternion multiplication is distributive over addition (), we can multiply the quaternions, "opening parenthesis". Note that quaternionic multiplication will be associative (), but in common case not commutative (there are such that ).

Quaternions can be also represented as , where — is vector of 3D space with unit vectors . Components and are called scalar and vector parts of quaternions. Assume we have quaternions and . Then their multiplication can be represented as . Consider multiplication of two quaternions with zero scalar parts.

.

Note that it can be shortly written as . Thus we have following formula for quaternionic multiplication:

Finally note that since , quaternion can be represented in form , where — are complex numbers. Also note that since , we have identity or, introducing notation : . Thus multiplication of quaternions and can be represented as . Number , introduced earlier, is called conjugated to complex number . So we have another one formula for quaternion multiplication:

Let's show that any non-zero quaternion is invertible, i.e. for each quaternion there is inverse such that . In this case we can divide quaternions by multiplicating with inverse element. Like in complex numers we can introduce for quaternion conjugated quaternion . From that formula we can see that . Thus we can introduce quaternionic norm (which is generalisation of length): and inverse element . Note that . Indeed from multiplication formula we get . From here we can see and .

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 ^[1] . This identity is equivalent to the following: . Let , then we can get . Considering scalar and vector parts separately we will get:

Note that because ||ab|| = ||a||·||b|| we can get that norms of vectors and are equal. First equation of the system means that orthogonal projections on axis are equal for and . Set of vectors with same norm and projection on some direction form a circumference around that direction. Thus we see that can be obtained from by its rotation around 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 . If it's not true, we can divide by square root of its norm, it will not affect . Now we can assume that , where . Also using first equation and using notation and , we can transform second equation in following way: . (On the other hand we subtracted from both sides, on the other one since we can see that )

Note that and are lying in the plane orthogonal to , thus, this is the rotation plane. Now we should see that is vector rotated around vector clockwise. Thus left part of second equation is , rotated around clockwise on angle , and the right one is , rotated around counter-clockwise on . Thus we see that is rotated around counter-clockwise on .

[1] . Scalar part of this product equals to .

And now the implementation. Ideone: #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;
}