Generalized N-Dimensional Rotors & Decompositions

[CaseySanchez]

Hello,

Given that rotors are expressed as shown below with their corresponding GA:
G^{2}: R = R_{0} * 1 + R_{12} * (e1 ^ e2)
G^{3}: R = R_{0} * 1 + R_{12} * (e1 ^ e2) + R_{13} * (e1 ^ e3) + R_{23} * (e2 ^ e3)
G^{4}: R = R_{0} * 1 + R_{12} * (e1 ^ e2) + R_{13} * (e1 ^ e3) + R_{14} * (e1 ^ e4) + R_{23} * (e2 ^ e3) + R_{24} * (e2 ^ e4) + R_{34} * (e3 ^ e4) + R_{1234} * (e1 ^ e2 ^ e3 ^ e4)

How can rotors be generalized into a geometric algebra of arbitrary dimensionality?

I would like to generalize these functions:

template<typename T, typename RotorExpression, typename TranslatorExpression>
auto decompose(clifford_expression<T, RotorExpression> const &rotor, clifford_expression<T, TranslatorExpression> const &translator)
{
    Vector<T> const rotor_vector = sp(rotor, static_cast<T>(-1.0) * (e1 ^ e2)) * e1 + sp(rotor, static_cast<T>(-1.0) * (e1 ^ e3)) * e2 + sp(rotor, static_cast<T>(-1.0) * (e2 ^ e3)) * e3;
    Vector<T> const rotor_project = project(rotor_vector, translator);

    Rotor<T> const twist = unit(sp(rotor, c<1>) + sp(rotor_project, e1) * (e1 ^ e2) + sp(rotor_project, e2) * (e1 ^ e3) + sp(rotor_project, e3) * (e2 ^ e3));
    Rotor<T> const swing = apply_rotor(rotor, inv(twist));
    
    return std::make_tuple(swing, twist);
}

template<typename T, typename RotorExpression>
auto constrain(T const &theta, clifford_expression<T, RotorExpression> const &rotor)
{
    T const norm = std::sin(theta * static_cast<T>(0.5));
    T const norm_sqr = norm * norm;

    Vector<T> const rotor_vector = sp(rotor, static_cast<T>(-1.0) * (e1 ^ e2)) * e1 + sp(rotor, static_cast<T>(-1.0) * (e1 ^ e3)) * e2 + sp(rotor, static_cast<T>(-1.0) * (e2 ^ e3)) * e3;

    T const rotor_vector_norm_sqr = rnorm_sqr(rotor_vector);

    if (rotor_vector_norm_sqr > norm_sqr) {
        Vector<T> const rotor_vector_normal = rotor_vector / std::sqrt(rotor_vector_norm_sqr);

        T const sign = sp(rotor, c<1>) < static_cast<T>(0.0) ? static_cast<T>(-1.0) : static_cast<T>(+1.0);
        
        return std::sqrt(1.0 - norm_sqr) * sign + sp(rotor_vector_normal, e1) * (e1 ^ e2) + sp(rotor_vector_normal, e2) * (e1 ^ e3) + sp(rotor_vector_normal, e3) * (e2 ^ e3);
    }
    
    return sp(rotor, c<1>) + sp(rotor, static_cast<T>(-1.0) * (e1 ^ e2)) * (e1 ^ e2) + sp(rotor, static_cast<T>(-1.0) * (e1 ^ e3)) * (e1 ^ e3) + sp(rotor, static_cast<T>(-1.0) * (e2 ^ e3)) * (e2 ^ e3);
}

Just as a quaternion can be decomposed into {3 \choose 2} - 1 = 2 quaternions about a vector, I presume an octonion can be decomposed into {4 \choose 2} - 1 = 5 octonions about a vector. Does this sound correct?

[CaseySanchez]

I have managed to find this paper specifically about swing-twist decomposition of quaternions:
https://arxiv.org/pdf/1506.05481.pdf

However, looking at the algorithm on p.22 it employs similar goofiness that had previously confused Dr. Fernandes.

[CaseySanchez]

This works if I manually extract the bivector of the rotor as is done in the comment. When I use checked_trivial_copy to extract the bivector component it throws ga::bad_checked_copy_exception:

DELETED

Main:

DELETED

I am still researching how to extract the remaining rotors in >3 dimensions.

[laffernandes]

To extract one step of the multivector you have to use take_grade(arg, k). The checked_trivial_copy procedure makes a copy, i.e., the target multivector must be capable of storing all non-zero coefficients of the source multivector.

[CaseySanchez]

Disregard that comment I made, it was a false-positive. I will correct and update it as soon as I can.

Thank you for the information regarding take_grade. I will use that instead.

[CaseySanchez]

I apologize for the previous false alarm Dr. Fernandes, I can confirm that the code below does in fact work as expected. I've tried to adhere to your GATL-isms, if you will:

template<typename CoefficientType, typename RotorExpression, typename VectorExpression>
auto decompose(clifford_expression<CoefficientType, RotorExpression> const &rotor, clifford_expression<CoefficientType, VectorExpression> const &vector)
{
    auto const rotor_bivector = take_grade(rotor, 2);
    auto const rotor_vector = dual(rotor_bivector);
    auto const rotor_project = project(rotor_vector, vector);

    auto const twist = unit(take_grade(rotor, 0) + undual(rotor_project));
    auto const swing = rotor * ~twist;
    
    return std::make_tuple(swing, twist);
}

A test confirming that swing * twist = rotor:

auto const rotor = 0.331029 + 0.849073 * (e1 ^ e2) - 0.309207 * (e1 ^ e3) + 0.27204 * (e2 ^ e3);
auto const vector = 1.0 * e1 + 1.0 * e2 + 1.0 * e3;

auto const [swing, twist] = decompose(rotor, vector);

std::cout << "rotor = " << rotor << std::endl;
std::cout << "vector = " << vector << std::endl;
std::cout << "swing = " << swing << std::endl; 
std::cout << "twist = " << twist << std::endl;
std::cout << "swing * twist = " << (swing * twist) << std::endl;

[laffernandes]

I apologize for the previous false alarm Dr. Fernandes,

Don't worry.

I can confirm that the code below does in fact work as expected. I've tried to adhere to your GATL-isms, if you will:

template<typename CoefficientType, typename RotorExpression, typename VectorExpression>
auto decompose(clifford_expression<CoefficientType, RotorExpression> const &rotor, clifford_expression<CoefficientType, VectorExpression> const &vector)
{
    auto const rotor_bivector = take_grade(rotor, 2);
    auto const rotor_vector = dual(rotor_bivector);
    auto const rotor_project = project(rotor_vector, vector);

    auto const twist = unit(take_grade(rotor, 0) + undual(rotor_project));
    auto const swing = rotor * ~twist;
    
    return std::make_tuple(swing, twist);
}

A test confirming that swing * twist = rotor:

auto const rotor = 0.331029 + 0.849073 * (e1 ^ e2) - 0.309207 * (e1 ^ e3) + 0.27204 * (e2 ^ e3);
auto const vector = 1.0 * e1 + 1.0 * e2 + 1.0 * e3;

auto const [swing, twist] = decompose(rotor, vector);

std::cout << "rotor = " << rotor << std::endl;
std::cout << "vector = " << vector << std::endl;
std::cout << "swing = " << swing << std::endl; 
std::cout << "twist = " << twist << std::endl;
std::cout << "swing * twist = " << (swing * twist) << std::endl;

Use take_grade(rotor, c<2>) instead of take_grade(rotor, 2). This way the compile code will be slightly faster.