Rotor Inaccuracies - Exp Bug

[CaseySanchez]

Hello Dr. Fernandes,

In the example below I rotate the "right" (+x), "forward" (+y), and "up" (+z) vectors by a 90-degree rotor about the e1^e2 bivector.

#include <gatl/ga3h.hpp>

using namespace ga3h;

int main(int argc, char **argv)
{
    auto const right = 1.0f * e1 + 0.0f * e2 + 0.0f * e3;
    auto const forward = 0.0f * e1 + 1.0f * e2 + 0.0f * e3;
    auto const up = 0.0f * e1 + 0.0f * e2 + 1.0f * e3;

    auto const rotor = exp<float>((-3.141592f * 0.5f * 0.5f) * (e1 ^ e2));

    auto const right_rotated = rotor * right * inv(rotor);
    auto const forward_rotated = rotor * forward * inv(rotor);
    auto const up_rotated = rotor * up * inv(rotor);

    std::cout << "Right: ";

    for (float const &value : right_rotated.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;

    std::cout << "Forward: ";

    for (float const &value : forward_rotated.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;
    
    std::cout << "Up: ";

    for (float const &value : up_rotated.values()) {
        std::cout << value << " ";
    }

    std::cout << std::endl;    

    return 0;
}

This outputs are the following:

Right: 0.236973 0.971516 0 0 
Forward: -0.971516 0.236973 0 0 
Up: 0 0 1 0 

The rotated "up" is as expected, resulting in the same vector. However, the rotated "right" is not, where the expected result would be the "forward" (+y) vector of 0.0f * e1 + 1.0f * e2 + 0.0f * e3. While the coefficient of e2 is certainly close to 1.0f, the coefficient of e1 is much larger than it should be. The result is similar for "forward".

Am I doing something wrong here? I am still in the process of learning Geometric Algebra, so perhaps I am misunderstanding something.

Additionally, in the above example how would I replace all of the auto keywords with their explicit types? I understand that "right", "forward", and "up" should be of type full_kvector_t<float, 3, 1>, but I am uncertain what the type of "rotor" and the rotated vectors should be.

Regards,
Casey

[CaseySanchez]

I believe I had made a mistake in using inv instead of the reverse, such as the following:
rotor * right * ~rotor

Making this change, however, yields even stranger results:

Right: 0.38315 1.5708 0 0 
Forward: -1.5708 0.38315 0 0 
Up: 0 0 1.61685 0 

EDIT:
It also appears that "normalizing" the rotor yields the same result in the original post:

auto const r = exp<float>((-3.141592f * 0.5f * 0.5f) * (e1 ^ e2));
auto const n = std::sqrt(std::abs(*((r * conjugate(r)).values().begin())));
auto const rotor = r * (1.0f / n);

With the result being the following:

Right: 0.236973 0.971516 0 0 
Forward: -0.971516 0.236973 0 0 
Up: 0 0 1 0 

[laffernandes]

Hi @Tannz0rz, I believe you have found a bug in the exp function. The result is correct when the rotor is set to:

float alpha = 3.141592f * 0.5f;
auto rotor = cos(alpha / 2) - sin(alpha / 2) * (e1 ^ e2);

Please, find below a cleaner and commented version of the source code:

#include <gatl/ga3h.hpp>

using namespace ga3h;

int main(int argc, char **argv)
{
    float alpha = 3.141592f * 0.5f;
    auto const right = 1.0f * e1 + 0.0f * e2 + 0.0f * e3; // The function `euclidean_vector(1.0f, 0.0f, 0.0f)` is a shortcut for defining vectors with homogeneous coordinates equals to zero. It is faster because it applies lazy evaluation.
    auto const forward = 0.0f * e1 + 1.0f * e2 + 0.0f * e3;
    auto const up = 0.0f * e1 + 0.0f * e2 + 1.0f * e3;

    auto rotor_exp = exp((-alpha / 2) * (e1 ^ e2));
    std::cout << "Rotor using exp(): " << rotor_exp << " (Reverse Norm: " << rnorm(rotor_exp) << ")" << std::endl; // This rotor is not correct (it seems a bug in the exp function).
    std::cout << std::endl;

    auto rotor = cos(alpha / 2) - sin(alpha / 2) * (e1 ^ e2);
    std::cout << "Rotor: " << rotor << " (Reverse Norm: " << rnorm(rotor) << ")" << std::endl; // This rotor is correct.
    std::cout << std::endl;

    auto const right_rotated = apply_rotor(rotor, right); // Equivalent to `rotor * right * inv(rotor)` and `rotor * right * ~rotor`, but faster because it explores lazy evaluation.
    auto const forward_rotated = apply_rotor(rotor, forward);
    auto const up_rotated = apply_rotor(rotor, up);

    std::cout << "Right: " << right_rotated << std::endl; // Use the << operator to print multivectors.
    std::cout << "Forward: " << forward_rotated << std::endl;
    std::cout << "Up: " << up_rotated << std::endl;

    return EXIT_SUCCESS;
}

I will investigate and fix the cause of the problem. But at the moment I'm swamped with work and so I won't be able to dedicate myself to it until next weekend.

Thank you for pointing out this issue.

About the auto keyword, replace it may be a little tricky since the lazy evaluation process is supposed to return the multivector type that is best for each particular expression. For instance, the types of x_ and y_ are different, even representing the same vector:

float alpha = 3.141592f * 0.5f;
auto rotor = cos(alpha / 2) - sin(alpha / 2) * (e1 ^ e3);
auto x = 10.0f * e1;
auto x_ = apply_rotor(rotor, x); // x_ = 3.8147e-06 * <e1> + 10 * <e3>
auto y = 10.0f * e1 + 0.0f * e2;
auto y_ = apply_rotor(rotor, x); // y_ = 3.8147e-06 * <e1> + 0 * <e2> + 10 * <e3>

Notice that x_ does not store a coefficient for e2.

[CaseySanchez]

Understood regarding your workload.

Is there a way to provide these resulting rotated vectors with the explicit type of full_multivector_t<float, 3>? I understand that it would contain redundant zero-valued k-vectors.

[laffernandes]

The option in that case is to declare the variable as full_vector_t<...>, full_kvector_t<...>, full_multivector_t<...>, etc. , and call trivial_copy or checked_trivial_copy explicitly:

#include <gatl/ga3h.hpp>

using namespace ga3h;

int main(int argc, char **argv)
{
    float alpha = 3.141592f * 0.5f;
    auto rotor = cos(alpha / 2) - sin(alpha / 2) * (e1 ^ e3);
    auto x = 10.0f * e1;

    full_vector_t<float, 4> r_;
    trivial_copy(apply_rotor(rotor, x), r_);
    
    std::cout << r_ << std::endl;

    return EXIT_SUCCESS;
}

[CaseySanchez]

Could the trivial_copy perhaps be fashioned into an assignment operator overload? Such that the resulting expression would be: full_vector_t<float, 4> r_ = apply_rotor(rotor, x);

[laffernandes]

Ok, I will do (almost) that with the exp bug fix. But I will call trivial_copy in release mode and checked_trivial_copy in debug mode. I believe it is safer that way.

For performance purposes, I strongly suggest avoiding these copies, as it is possible to extract the coefficients (even the zeros) from the result inferred by auto.

[CaseySanchez]

One last semi-related question Dr. Fernandes. How would I approach the explicit type declaration of a rotor/quaternion, which is equivalent to declaring a multivector consisting of a scalar part and three bivector parts in G^{3}?

I understand your feelings with regard to the use of auto with your library in terms of best practices, I'm just trying to have a better understanding of the underlying types of the library itself.

EDIT: Furthermore, I would like to pass these rotor/quaternion types (for which I believe there are 8 possible permutations of multivectors in G^{3}, including the identity) to a member function using SFINAE e.g. std::enable_if_t<is_rotor<T>::value>. Do you have an example of such type constraints?

[laffernandes]

The dirty trick uses the decltype keyword:

using rotor_t = decltype(float() + float() * (e1^e2) + float() * (e1^e3) + float() * (e2^e3));

You can do that to deduce the resulting type of any C++ expression. It is part of the language, not the library. So, by defining an expression with the components you want, you get the multivector type.

The order of the elements is important for GATL. Using the expression presented above, you will get the only permutation that GATL accepts as valid. That is part of the library, so it can not be changed. It means that rotor1_t and rotor2_t below will be the same type. The convention of the library uses the order of the components declared in rotor1_t.

using rotor1_t = decltype(float() + float() * (e1^e2) + float() * (e1^e3) + float() * (e2^e3));
using rotor2_t = decltype(float() + float() * (e1^e2) + float() * (e2^e3) + float() * (e1^e3));

Using the std::enable_if_t<is_rotor<T>::value> test in compilation time will not work because a multivector consists of a scalar part plus three bivector parts in G^{3} may not be a rotor. This can only be checked at runtime.

[CaseySanchez]

I wanted to update you with what I have been doing exactly. I made a Transform class that allows for transform hierarchies using GATL, as seen in the Gist below:
https://gist.github.com/Tannz0rz/cbd2938e671990092e800ad8c57ba7a2

I deliberately keep the rotors and translators separate in case the user may wish to set a global rotor while retaining a local translator, for example.

As you can see there are a few instances where I do not believe I can avoid using checked_trivial_copy, where the previously-mentioned operator= overload would help. Moreover, a copy constructor employing this behavior would help too. When you do implement this as you said you would, to inform users of this you may wish to issue a #pragma message "checked_trivial_copy/trivial_copy is being used", or something similar.

All-in-all, I want to say thank you so much for this library Dr. Fernandes, I am really enjoying it so far. This is how math in C++ should be. It is time to move away from Eigen and onto GATL!

Could you perhaps review my code and make suggestions? I wish I could move the setParent code into the constructor, but that is not possible using shared_from_this. I would like to share this example on the Bivector.net Discord and on the r/cpp Subreddit and spread the word about GATL. Correct me if I'm wrong, but given your design philosophy I would argue that it is the most performant of all GA libraries currently in existence.