The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings Garamon
Summary Activity SCM

SCM Repository

authorbreuils <stephane.breuils@u-pem.frcd>
Tue, 30 Apr 2019 19:18:14 +0000 (21:18 +0200)
committerbreuils <stephane.breuils@u-pem.frcd>
Tue, 30 Apr 2019 19:18:14 +0000 (21:18 +0200)
README.md
data/CMakeLists.txt
data/Mvec.hpp
data/sample/src/main.cpp
src/Directory.cpp
src/ProductTools.hpp

index 2bf3acb..69baa44 100644 (file)
--- a/README.md
+++ b/README.md
@@ -1,68 +1,80 @@
-Garamon Generator
-=================
-
-Garamon (Geometric Algebra Recursive and Adaptative Monster) is a generator of C++ libraries dedicated to Geometric Algebra.
-From a configuration file, the library generator generate the source code of the specified Algebra as well as a compatible sample code and some documentation.
-
-## Features
-    * MIT Licence allows free use in all software (including GPL and commercial)
-    * multi-platform (Windows, Linux, Unix, Mac)
-
-
-Install
-=======
-
-## Dependencies
-    * 'Eigen 3.3.4'  or more [Eigen](http://eigen.tuxfamily.org)
-    * 'Cmake 3.10' or more
-    * (optional) 'Doxygen' if you want documentation pages
-
-## Compilator tested
-    * gcc 5.4.0
-    * clang 4
-    * apple-clang 900.0.39.2
-    * MinGW 7.2.0
-
-## Install for Linux-Mac from the terminal
-    * Check the dependencies
-    * From Garamon Generator root directory
-        * 'mkdir build'
-        * 'cd build'
-        * 'cmake ..'
-        * check that the cmake output has no errors
-        * 'make'
-        * the binary executable is on the 'build' directory
-        * (optional and not required) 'sudo make install'
-
-## Install for Windows with MinGW, using Windows Power Shell
-    * Check the dependencies
-    * From Garamon Generator root directory
-        * 'mkdir build'
-        * 'cd build'
-        * 'cmake -G "MinGW Makefiles" ..'
-        * check that the cmake output has no errors
-               * 'mingw32-make'
-        * the binary executable is on the 'build' directory
-        * (optional and not required) 'mingw32-make install'
-
-
-Usage
-=====
-
-## Generate a library
-
-    * define the algebra to generate: chose a configuration file (.conf) on the 'conf' directory or create your own.
-    * run the binary executable (from the 'build' directory) with the configuration file as argument
-      > ./garamon_generator file.conf
-    * the generated library is located in 'build/output' directory
-
-
-Notes
-=====
-
-## Authors
-    * Vincent Nozick (Universite Paris-Est Marne-la-Vallee, France)
-    * Stephane Breuils (Universite Paris-Est Marne-la-Vallee, France)
-
-## Contact
-    * vincent.nozick (at) u-pem.fr
+Garamon Generator\r
+=================\r
+\r
+Garamon (Geometric Algebra Recursive and Adaptative Monster) is a generator of C++ libraries dedicated to Geometric Algebra.\r
+From a configuration file, the library generator generate the source code of the specified Algebra as well as a compatible sample code and some documentation.\r
+\r
+## Features\r
+    * MIT Licence allows free use in all software (including GPL and commercial)\r
+    * multi-platform (Windows, Linux, Unix, Mac)\r
+\r
+\r
+Install\r
+=======\r
+\r
+## Dependencies\r
+    * 'Eigen 3.3.4'  or more [Eigen](http://eigen.tuxfamily.org)\r
+    * 'Cmake 3.10' or more\r
+    * (optional) 'Doxygen' if you want documentation pages\r
+\r
+## Compilator tested\r
+    * gcc 5.4.0\r
+    * clang 4\r
+    * apple-clang 900.0.39.2\r
+    * MinGW 7.2.0\r
+       * MSVC 19.14.26430.0 (Visual Studio 15.7.3)\r
+\r
+## Install for Linux-Mac from the terminal\r
+    * Check the dependencies\r
+    * From Garamon Generator root directory\r
+        * 'mkdir build'\r
+        * 'cd build'\r
+        * 'cmake ..'\r
+        * check that the cmake output has no errors\r
+        * 'make'\r
+        * the binary executable is on the 'build' directory\r
+        * (optional and not required) 'sudo make install'\r
+\r
+## Install for Windows with MinGW, using Windows Power Shell\r
+    * Check the dependencies\r
+    * From Garamon Generator root directory\r
+        * 'mkdir build'\r
+        * 'cd build'\r
+        * 'cmake -G "MinGW Makefiles" ..'\r
+        * check that the cmake output has no errors\r
+               * 'mingw32-make'\r
+        * the binary executable is on the 'build' directory\r
+        * (optional and not required) 'mingw32-make install'\r
+\r
+## Install for Windows with Visual Studio 15 2017 Win64, using Windows Power Shell or cmd\r
+    * Check the dependencies\r
+    * From Garamon Generator root directory\r
+        * 'mkdir build'\r
+        * 'cd build'\r
+        * 'cmake -G "Visual Studio 15 2017 Win64" ..'\r
+        * check that the cmake output has no errors\r
+               * open the file garamon_generator.sln with Visual Studio\r
+               * generate the project "ALL_BUILD" with Release configuration\r
+        * the binary executable is on the 'Release' directory\r
+\r
+\r
+Usage\r
+=====\r
+\r
+## Generate a library\r
+\r
+    * define the algebra to generate: chose a configuration file (.conf) on the 'conf' directory or create your own.\r
+    * run the binary executable (from the 'build' directory) with the configuration file as argument\r
+      > ./garamon_generator file.conf\r
+    * the generated library is located in 'build/output' directory\r
+\r
+\r
+Notes\r
+=====\r
+\r
+## Authors\r
+    * Vincent Nozick (Universite Paris-Est Marne-la-Vallee, France)\r
+    * Stephane Breuils (Universite Paris-Est Marne-la-Vallee, France)\r
+\r
+## Contact\r
+    * vincent.nozick (at) u-pem.fr\r
index 4055ef6..431553b 100644 (file)
@@ -1,54 +1,61 @@
-cmake_minimum_required(VERSION 3.5)
-
-project(cmake_project_name_original_case)
-set(CMAKE_BUILD_TYPE Release)
-
-# include eigen
-find_package(Eigen3 REQUIRED)
-if(${EIGEN3_FOUND})
-    message(STATUS "lib EIGEN3 found")
-    message(STATUS "  version " ${EIGEN3_VERSION_STRING})
-    message(STATUS "  include " ${EIGEN3_INCLUDE_DIR})
-else()
-    include_directories("/usr/include/eigen3") # manually specify the include location
-endif()
-
-# call the CMakeLists.txt to make the documentation (Doxygen)
-# > 'make html' to generate the documentation
-add_subdirectory(doc)
-
-# compilation flags
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -O2 -std=c++14 -Wno-return-local-addr")
-set(CMAKE_CXX_FLAGS_DEBUG_INIT "-g")
-#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ftime-report ") # comilation profiling:  -Q -ftime-report
-#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fno-elide-constructors") # force the compilor to use our move 'copy constructor'
-
-# files to compile
-file(GLOB_RECURSE source_files src/cmake_project_name_original_case/*.cpp)
-file(GLOB_RECURSE header_files src/cmake_project_name_original_case/*.hpp)
-
-# display info
-message(STATUS "  sources")
-foreach(src_file ${source_files})
-    message("          " ${src_file})
-endforeach()
-message(STATUS "  headers")
-foreach(src_file ${header_files})
-    message("          " ${src_file})
-endforeach()
-
-# generate a library
-add_library(cmake_project_name_original_case SHARED ${source_files} ${header_files})
-#add_library(cmake_project_name_original_case STATIC ${source_files} ${header_files})
-
-# include directory path
-include_directories(src)
-include_directories(${EIGEN3_INCLUDE_DIR})
-
-# install lib
-install(FILES ${header_files} ${source_files} DESTINATION /usr/local/include/cmake_project_name_original_case)
-install(TARGETS cmake_project_name_original_case
-        RUNTIME DESTINATION /usr/local/bin
-        LIBRARY DESTINATION /usr/local/lib
-        ARCHIVE DESTINATION /usr/local/lib)
-
+cmake_minimum_required(VERSION 3.5)\r
+\r
+project(cmake_project_name_original_case)\r
+set(CMAKE_BUILD_TYPE Release)\r
+\r
+# include eigen\r
+find_package(Eigen3 REQUIRED)\r
+if(${EIGEN3_FOUND})\r
+    message(STATUS "lib EIGEN3 found")\r
+    message(STATUS "  version " ${EIGEN3_VERSION_STRING})\r
+    message(STATUS "  include " ${EIGEN3_INCLUDE_DIR})\r
+else()\r
+    include_directories("/usr/include/eigen3") # manually specify the include location\r
+endif()\r
+\r
+# call the CMakeLists.txt to make the documentation (Doxygen)\r
+# > 'make html' to generate the documentation\r
+add_subdirectory(doc)\r
+\r
+# compilation flags\r
+if (MSVC)\r
+       set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -O2 -std=c++14")\r
+else()\r
+       set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -O2 -std=c++14 -Wno-return-local-addr")\r
+endif()\r
+set(CMAKE_CXX_FLAGS_DEBUG_INIT "-g")\r
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ftime-report ") # comilation profiling:  -Q -ftime-report\r
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fno-elide-constructors") # force the compilor to use our move 'copy constructor'\r
+\r
+# files to compile\r
+file(GLOB_RECURSE source_files src/cmake_project_name_original_case/*.cpp)\r
+file(GLOB_RECURSE header_files src/cmake_project_name_original_case/*.hpp)\r
+\r
+# display info\r
+message(STATUS "  sources")\r
+foreach(src_file ${source_files})\r
+    message("          " ${src_file})\r
+endforeach()\r
+message(STATUS "  headers")\r
+foreach(src_file ${header_files})\r
+    message("          " ${src_file})\r
+endforeach()\r
+\r
+# generate a library\r
+if (MSVC)\r
+       add_library(cmake_project_name_original_case STATIC ${source_files} ${header_files})\r
+else()\r
+       add_library(cmake_project_name_original_case SHARED ${source_files} ${header_files})\r
+endif()\r
+\r
+# include directory path\r
+include_directories(src)\r
+include_directories(${EIGEN3_INCLUDE_DIR})\r
+\r
+# install lib\r
+install(FILES ${header_files} ${source_files} DESTINATION /usr/local/include/cmake_project_name_original_case)\r
+install(TARGETS cmake_project_name_original_case\r
+        RUNTIME DESTINATION /usr/local/bin\r
+        LIBRARY DESTINATION /usr/local/lib\r
+        ARCHIVE DESTINATION /usr/local/lib)\r
+\r
index a466fe5..78485d6 100644 (file)
-// Copyright (c) 2018 by University Paris-Est Marne-la-Vallee
-// Mvec.hpp
-// This file is part of the Garamon for project_namespace.
-// Authors: Stephane Breuils and Vincent Nozick
-// Contact: vincent.nozick@u-pem.fr
-//
-// Licence MIT
-// A a copy of the MIT License is given along with this program
-
-/// \file Mvec.hpp
-/// \author Stephane Breuils, Vincent Nozick
-/// \brief Class to define a multivector and its basic operators in the Geometric algebra of project_namespace.
-
-
-// Anti-doublon
-#ifndef project_inclusion_guard
-#define project_inclusion_guard
-#pragma once
-
-// External Includes
-#include <Eigen/Core>
-#include <list>
-#include <iostream>
-#include <cmath>
-#include <limits>
-
-// Internal Includes
-#include "project_namespace/Utility.hpp"
-#include "project_namespace/Constants.hpp"
-
-#include "project_namespace/Outer.hpp"
-#include "project_namespace/Inner.hpp"
-#include "project_namespace/Geometric.hpp"
-
-#include "project_namespace/OuterExplicit.hpp"
-#include "project_namespace/InnerExplicit.hpp"
-#include "project_namespace/GeometricExplicit.hpp"
-
-/*!
- * @namespace project_namespace
- */
-namespace project_namespace{
-
-
-    /// \class Kvec
-    /// \brief class defining a single grade component of a multivector.
-    template<class T>
-    struct Kvec{
-
-        Eigen::Matrix<T, Eigen::Dynamic, 1> vec;  /*!< dynamic vector of Eigen Library */
-
-        unsigned int grade; /*!< grade k of the k-vector */
-
-        /// \brief operator == to test the equality between 2 k-vectors.
-        bool operator == (const Kvec<T>& other) const {
-            return vec == other.vec;
-        }
-    };
-
-
-
-    /// \class Mvec
-    /// \brief class defining multivectors.
-    template<typename T = double>
-    class Mvec {
-
-    protected:
-        std::list<Kvec<T>> mvData;  /*!< set of k-vectors, mapped by grade */
-        unsigned int gradeBitmap;   /*!< ith bit to 1 if grade i is contained in the multivector */
-
-    public:
-
-        /// \brief Default constructor, generate an empty multivector equivalent to the scalar 0.
-        Mvec();
-
-        /// \brief Copy constructor
-        /// \param mv - the multivector to be copied
-        Mvec(const Mvec& mv);
-
-        /// \brief Copy constructor of Mvec
-        /// \param mv - the multivector which has to be copied
-        Mvec(Mvec&& mv); // move constructor
-
-
-        template <typename U>
-        friend class Mvec;
-
-        /// \brief Copy constructor of Mvec from different types
-        /// \param mv - the multivector with another template type
-        template<typename U>
-        Mvec<T>(const Mvec<U> &mv);
-
-        /// \brief Constructor of Mvec from a scalar
-        /// \param val - scalar value
-        template<typename S>
-        Mvec(const S val);
-
-        /// Destructor
-        ~Mvec();
-
-        /// \brief Overload the assignment operator. No need any copy when the argument is an R-value, just need to move.
-        /// \param mv - Mvec used as a Rvalue
-        /// \return assign other to this
-        Mvec& operator=(Mvec&& mv);
-
-        /// \brief Overload the assignment operator
-        /// \param mv - Mvec
-        /// \return assign mv to this object
-        Mvec& operator=(const Mvec& mv);
-
-        /// \brief defines the addition between two Mvec
-        /// \param mv2 - second operand of type Mvec
-        /// \return this + mv2
-        Mvec operator+(const Mvec &mv2) const;
-
-        /// \brief defines the addition between a Mvec and a scalar
-        /// \param value - second operand (scalar)
-        /// \return this + scalar
-        template<typename S>
-        Mvec operator+(const S &value) const;
-
-        /// \brief defines the addition between a scalar and a Mvec
-        /// \param value a scalar
-        /// \param mv - second operand of type Mvec
-        /// \return scalar + mv
-        template<typename U, typename S>
-        friend Mvec<U> operator+(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the += operator, corresponds to this += mv
-        /// \param mv - Mvec to be added to this object
-        /// \return this += mv
-        Mvec& operator+=(const Mvec& mv);
-
-        /// \brief defines the opposite of a multivector
-        /// \param mv: operand of type Mvec
-        /// \return -mv
-        template<typename U>
-        friend Mvec<U> operator-(const Mvec<U> &mv); // unary operator -mv1
-
-        /// \brief defines the difference between two Mvec
-        /// \param mv2 - second operand of type Mvec
-        /// \return this - mv2
-        Mvec<T> operator-(const Mvec<T> &mv2) const;
-
-        /// \brief defines the difference between a Mvec and a scalar
-        /// \param value - second operand (scalar)
-        /// \return this - scalar
-        template<typename S>
-        Mvec operator-(const S &value) const;
-
-        /// \brief defines the difference between a scalar and a Mvec
-        /// \param value a scalar
-        /// \param mv - second operand of type Mvec
-        /// \return scalar - mvX
-        template<typename U, typename S>
-        friend Mvec<U> operator-(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the -= operator, corresponds to this -= mv
-        /// \param mv - Mvec to be added to this object
-        /// \return this -= mv
-        Mvec& operator-=(const Mvec& mv);
-
-        /// \brief defines the outer product between two multivectors
-        /// \param mv2 - a multivector
-        /// \return this^mv2
-        Mvec operator^(const Mvec &mv2) const;
-
-        /// \brief defines the outer product a multivector and a scalar
-        /// \param value - a scalar
-        /// \return this^value
-        template<typename S>
-        Mvec operator^(const S &value) const;
-
-        /// \brief defines the outer product between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return value ^ mv
-        template<typename U, typename S>
-        friend Mvec<U> operator^(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the outer product with operator =, corresponds to this ^= mv
-        /// \param mv - Mvec to be wedged to this object
-        /// \return this ^= mv
-        Mvec& operator^=(const Mvec& mv);
-
-        /// \brief defines the inner product between two multivectors
-        /// \param mv2 - a multivector
-        /// \return this.mv2
-        Mvec operator|(const Mvec &mv2) const;
-
-        /// \brief defines the inner product between a multivector and a scalar
-        /// \param value - a scalar
-        /// \return this.value = 0 (empty multivector)
-        template<typename S>
-        Mvec operator|(const S &value) const;
-
-        /// \brief defines the inner product between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return mv.value = 0 (empty multivector)
-        template<typename U, typename S>
-        friend Mvec<U> operator|(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the inner product with operator =, corresponds to this |= mv
-        /// \param mv - Mvec to be computed in the inner product with this object
-        /// \return this |= mv
-        Mvec& operator|=(const Mvec& mv);
-
-        /// \brief defines the right contraction between two multivectors
-        /// \param mv2 - a multivector
-        /// \return the right contraction : $this \\lfloor mv2$
-        Mvec operator>(const Mvec &mv2) const;
-
-        /// \brief defines the right contraction between a multivector and a scalar
-        /// \param value - a scalar
-        /// \return $this \\lfloor value$
-        template<typename S>
-        Mvec operator>(const S &value) const;
-
-        /// \brief defines the right contraction between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return $value \\lfloor this = 0$ (empty multivector)
-        template<typename U, typename S>
-        friend Mvec<U> operator>(const S &value, const Mvec<U> &mv);
-
-        /// \brief defines the left contraction between two multivectors
-        /// \param mv2 - a multivector
-        /// \return the left contraction $this \\rfloor mv2$
-        Mvec operator<(const Mvec &mv2) const;
-
-        /// \brief defines the left contraction between a multivector and a scalar
-        /// \param value - a scalar
-        /// \return $value \\rfloor this$ = 0 (empty multivector)
-        template<typename S>
-        Mvec operator<(const S &value) const;
-
-        /// \brief defines the left contraction between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return $this \\rfloor value$
-        template<typename U, typename S>
-        friend Mvec<U> operator<(const S &value, const Mvec<U> &mv);
-
-        /// \brief defines the geometric product between two multivectors
-        /// \param mv2 - a multivector
-        /// \return mv1*mv2
-        Mvec operator*(const Mvec &mv2) const;
-
-    project_singular_metric_comment_begin
-        /// \brief defines the outer product between two multivectors, where the second multivector is dualized during the product computation.
-        /// \param mv2 - a multivector that will be dualized during the wedge with the calling multivector.
-        /// \return a multivector.
-        Mvec<T> outerPrimalDual(const Mvec<T> &mv2) const;
-
-        /// \brief defines the outer product between two multivectors, where the first multivector is dualized during the product computation.
-        /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.
-        /// \return a multivector.
-        Mvec<T> outerDualPrimal(const Mvec<T> &mv2) const;
-
-        /// \brief defines the outer product between two multivectors, where both multivectors are dualized during the product computation.
-        /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.
-        /// \return a multivector.
-        Mvec<T> outerDualDual(const Mvec<T> &mv2) const;
-
-    project_singular_metric_comment_end
-        
-
-        /// \brief defines the scalar product between two multivectors (sum of the inner products between same grade pairs from the 2 multivectors)
-        /// \param mv2 - a multivector
-        /// \return a scalar
-        Mvec<T> scalarProduct(const Mvec<T> &mv2) const;
-
-        /// \brief defines the Hestenes product between two multivectors (Inner product - scalar product)
-        /// \param mv2 - a multivector
-        /// \return a multivector.
-        Mvec<T> hestenesProduct(const Mvec<T> &mv2) const;
-
-        /// \brief defines the geometric product between a multivector and a scalar
-        /// \param value - a scalar
-        /// \return mv2*value
-        template<typename S>
-        Mvec operator*(const S &value) const;
-
-        /// \brief defines the geometric product between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return value*mv2
-        template<typename U, typename S>
-        friend Mvec<U> operator*(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the geometric product with operator =, corresponds to this *= mv
-        /// \param mv - Mvec to be multiplied to this object
-        /// \return this *= mv
-        Mvec& operator*=(const Mvec& mv);
-
-        /// \brief defines the geometric product with a multivector and the inverse of a second multivector
-        /// \param mv2 - a multivector
-        /// \return this / mv2
-        Mvec operator/(const Mvec &mv2) const;
-
-        /// \brief defines the scalar inverse between a multivector and a scalar
-        /// \param value - a scalar
-        /// \return this / value
-        template<typename S>
-        Mvec operator/(const S &value) const;
-
-        /// \brief defines the inverse product between a scalar and a multivector
-        /// \param value - a scalar
-        /// \param mv - a multivector
-        /// \return value / mv
-        template<typename U, typename S>
-        friend Mvec<U> operator/(const S &value, const Mvec<U> &mv);
-
-        /// \brief Overload the inverse with operator =, corresponds to this /= mv
-        /// \param mv - Mvec to be inversed to this object
-        /// \return this /= mv
-        Mvec& operator/=(const Mvec& mv);
-
-        /// \brief the reverse of a multivector, i.e. if mv = a1^a2^...^an, then reverse(mv) = an^...^a2^a1
-        /// \param mv - a multivector
-        /// \return reverse of mv
-        template<typename U>
-        friend Mvec<U> operator~(const Mvec<U> &mv);
-
-    project_singular_metric_comment_begin
-        /// \brief the dual of a k-vector is defined as $A_k^* = A_k \\lcont I_n^{-1}$, for a multivector, we just dualize all its components.
-        /// \param mv - a multivector
-        /// \return the dual of mv
-        template<typename U>
-        friend Mvec<U> operator!(const Mvec<U> &mv);
-    project_singular_metric_comment_end
-        /// \brief boolean operator that tests the equality between two Mvec
-        /// \param mv2 - second operand of type Mvec
-        /// \return whether two Mvec have the same coefficients
-        inline bool operator==(const Mvec& mv2){
-            if(gradeBitmap != mv2.gradeBitmap)
-                return false;
-
-            return mvData == mv2.mvData;   //// listcaca : ca marche que si les listes sont ordonnees
-        }
-
-        /// \brief compute the inverse of a multivector
-        /// \return - the inverse of the current multivector
-        Mvec<T> inv() const;
-
-
-            /// \brief operator to test whether two Mvec have not the same coefficients
-        /// \param mv2 - second operand of type Mvec
-        /// \return boolean that specify the non-equality between two Mvec
-        inline bool operator!=(const Mvec& mv2){
-            return !(mvData == mv2.mvData);
-        }
-
-        /// \brief Display all the non-null basis blades of this objects
-        /// \param stream - destination stream
-        /// \param mvec - the multivector to be outputed
-        /// \return a stream that contains the list of the non-zero element of the multivector
-        template<typename U>
-        friend std::ostream& operator<< (std::ostream& stream, const Mvec<U> &mvec);
-
-        /// \brief overload the casting operator, using this function, it is now possible to compute : float a = float(mv);
-        /// \return the scalar part of the multivector
-        operator T () {
-            // if no scalar part, return
-
-            if( (gradeBitmap & 1) == 0 )
-                return 0;
-
-            // assuming now that the scalar part exists, return it
-            return findGrade(0)->vec[0];
-        }
-
-        /// \brief Overload the [] operator to assign a basis blade to a multivector. As an example, float a = mv[E12] = 42.
-        /// \param idx - the basis vector index related to the query.
-        /// \return the coefficient of the multivector corresponding to the "idx" component.
-        /// \todo handle the cases when the number of indices is higher than the dimension, and when the index is too high
-        T& operator[](const int idx ){
-
-            const unsigned int grade = xorIndexToGrade[idx];
-            const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];
-
-            auto it = mvData.begin();
-            while(it != mvData.end()){
-
-                // if grade not reach yet, continue
-                if(it->grade < grade)
-                    ++it;
-                else{
-
-                    // if grade found
-                    if(it->grade == grade)
-                        return it->vec[idxHomogeneous];
-
-                    // if grade exceed, create it and inster it before the current element
-                    Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};
-                    auto it2 = mvData.insert(it,kvec);
-                    gradeBitmap |= 1 << (grade);
-                    return it2->vec[idxHomogeneous];
-                }
-            }
-
-            // if the searched element should be added at the end, add it
-            Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};
-
-            auto it2 = mvData.insert(mvData.end(),kvec);
-            gradeBitmap |= 1 << (grade);
-            return it2->vec[idxHomogeneous];
-        }
-
-        /// \brief Overload the [] operator to copy a basis blade of this multivector. As an example, float a = mv[E12].
-        /// \param idx - the basis vector index related to the query.
-        /// \return the coefficient of the multivector corresponding to the "idx" component.
-        const T& operator[](const int idx) const{
-
-            const unsigned int grade = xorIndexToGrade[idx];
-            const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];
-
-            auto it = mvData.begin();
-            while(it != mvData.end()){
-
-                // if grade not reach yet, continue
-                if(it->grade < grade)
-                    ++it;
-                else{
-
-                    // if grade found
-                    if(it->grade == grade)
-                        return it->vec[idxHomogeneous];
-
-                    // if grade exceed, return a reference on ... humm ...
-                    T val = T(0);
-                    return val;
-                }
-            }
-
-            // if the searched element should be added at the end, humm ...
-            T val = T(0);   /// \todo find an elegant way to correct this
-            return val;
-        }
-
-/*
-        /// \cond DEV
-        // Variadic operators
-        /// \brief Overload the () operator to assign a basis blade to a multivector. As an example, consider we have a Mvec a and we want to set its e23 component to 4.7, using operator(), it will consists in writing a[2,3]=4.8.
-        /// \tparam Arguments - denotes a variadic list of index
-        /// \param listIndices - the list of indices
-        /// \return the right index in the homogeneous vectorXd
-        /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high
-        template<typename... List>
-        T& operator()(List ...listIndices){
-            // operation on these arguments
-            // First determine the grade, simply use the variadic function sizeof...()
-            // This function is correct in terms of indexing
-            const int gd = sizeof...(listIndices); //
-
-            // Second, from these parameters, compute the index in the corresponding VectorXd
-            const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);//Binomial<algebraDimension,binomCoef>::value;//binomCoef - sumB(first,listIndices...);//computeOrderingIndex<>(listIndices...);//(Binomial<algebraDimension,gd>::value-1);//-computeOrderingIndex<int,List>(4,2,listIndices...);//computeOrderingIndex<algebraDimension,gd>(listIndices...);//idxVariadic<algebraDimension,gd,first,listIndices...>::value; //computeOrderingIndex<algebraDimension,gd,listIndices...>();//idxVariadic<algebraDimension,gd,listIndices...>::value; // (Binomial<algebraDimension,gd>::value-1)- idxVariadic<algebraDimension,gd,listIndices...>::value
-
-            if(mvData.count(gd)>0)
-                return (mvData.at(gd))(idx);
-
-            mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);
-
-            return mvData[gd](idx);
-        }
-        /// \endcond
-*/
-
-        /// \cond DEV
-        /// \brief returns a multivector with one component to 1
-        /// \param grade : grade of the component to enable
-        /// \param index : index of the parameter in the k-vector (k = grade)
-        inline Mvec componentToOne(const unsigned int grade, const int index){
-            Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
-            kvec.vec[index] = T(1);
-            Mvec mv1;
-            mv1.mvData.push_back(kvec);
-            mv1.gradeBitmap = 1 << (grade);
-
-            return mv1;
-        }
-        /// \endcond // do not comment this functions
-
-        /// \cond DEV
-        /// \brief create a VectorXd if it has not yet been created
-        /// \param grade - grade of the considered kvector
-        /// \return nothing
-        inline typename std::list<Kvec<T>>::iterator createVectorXdIfDoesNotExist(const unsigned int grade){
-
-            auto it = mvData.begin();
-            while(it != mvData.end()){
-
-                // if grade not reach yet, continue
-                if(it->grade < grade)
-                    ++it;
-                else{
-
-                    // if grade found
-                    if(it->grade == grade)
-                        return it;
-
-                    // if grade exceed, create it and inster it before the current element
-                    Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
-                    auto it2 = mvData.insert(it,kvec);
-                    gradeBitmap |= 1 << (grade);
-                    return it2;
-                }
-            }
-
-            // if the searched element should be added at the end, add it
-            Kvec<T> kvec = {.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),.grade=grade};
-            auto it2 = mvData.insert(mvData.end(),kvec);
-            gradeBitmap |= 1 << (grade);
-            return it2;
-        }
-
-
-        /// \cond DEV
-        /// \brief modify the element of the multivector whose grade is "grade"
-        /// \param grade : the grade to access
-        /// \param indexVectorXd : index of the k-vector (with k = "grade")
-        /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd
-        inline T& at(const int grade, const int indexVectorXd){
-            auto it = createVectorXdIfDoesNotExist(grade);
-            return it->vec.coeffRef(indexVectorXd);  
-        }
-        /// \endcond // do not comment this functions
-
-        /// \cond DEV
-        /// \brief return the element of the multivector whose grade is "grade"
-        /// \param grade : the grade to access
-        /// \param indexVectorXd : index of the k-vector (with k = "grade")
-        /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd
-        inline T at(const int grade, const int indexVectorXd) const{
-
-            if((gradeBitmap & (1<<(grade))) == 0 )
-                return T(0);
-
-            auto it = findGrade(grade);
-            return it->vec.coeff(indexVectorXd);
-        }
-        /// \endcond // do not comment this functions
-
-        /// \brief the L2-norm of the mv is sqrt( abs( mv.mv ) )
-        /// \return the L2-norm of the multivector (as a double)
-        T inline norm() const {
-            return sqrt( fabs( (*this).scalarProduct( this->reverse() ) ));
-        };
-
-        /// \brief the L2-norm over 2 of the mv is mv.mv
-        /// \return the L2-norm of the multivector (as a double)
-        T inline quadraticNorm() const {
-            return ( this->reverse() | (*this) );
-        };
-
-    project_singular_metric_comment_begin
-        /// \brief compute the dual of a multivector (i.e mv* = reverse(mv) * Iinv)
-        /// \return - the dual of the multivector
-        Mvec<T> dual() const;
-    project_singular_metric_comment_end
-        /// \brief compute the reverse of a multivector
-        /// \return - the reverse of the multivector
-        Mvec<T> reverse() const;
-
-        /// \brief search in the multivector for a Kvec of grade "grade"
-        /// \return return a const iterator on the Kvec if exist, else return mvData.end()
-        inline typename std::list<Kvec<T>>::const_iterator findGrade(const unsigned int & gradeToFind) const {
-            for(auto it = mvData.begin();it != mvData.end();++it){
-                if(it->grade==gradeToFind){
-                    return it;
-                }
-            }
-            return mvData.end();
-        }
-
-        /// \brief search in the multivector for a Kvec of grade "grade"
-        /// \return return an iterator on the Kvec if exist, else return mvData.end()
-        inline typename std::list<Kvec<T>>::iterator findGrade(const unsigned int & gradeToFind) {
-            for(auto it = mvData.begin();it != mvData.end();++it){
-                if(it->grade==gradeToFind){
-                    return it;
-                }
-            }
-            return mvData.end();
-        }
-
-        /// \brief return the (highest) grade of the multivector
-        /// \return the highest grade of the multivector
-        inline int grade() const {
-            return (mvData.begin() == mvData.end()) ? 0 : mvData.rbegin()->grade;
-        }
-
-        /// \brief return the all non-zero grades of the multivector (several grades for non-homogeneous multivectors)
-        /// \return a list of the grades encountered in the multivector
-        std::vector<unsigned int> grades() const;
-
-        /// \brief returns a multivector that contains all the components of this multivector whose grade is i
-        /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.
-        Mvec grade(const int i) const;
-
-        /// \brief tell whether the multivector has grade component
-        /// \param grade - grade of the considered kvector
-        /// \return true if multivector has grade component, false else
-        inline bool isGrade(const unsigned int grade) const{
-            return ( (gradeBitmap & (1<<(grade+1))) != 0 );
-        }
-
-        /// \brief partially or completely erase the content of a multivector
-        /// \param grade : if < 0, erase all the multivector, else just erase the part of grade "grade".
-        void clear(const int grade = -1);
-
-        /// \brief check is a mutivector is empty, i.e. corresponds to 0.
-        /// \return True if the multivector is empty, else False.
-        inline bool isEmpty() const {
-            return mvData.empty();
-        }
-
-        /// \brief A multivector is homogeneous if all its components have the same grade.
-        /// \return True if the multivector is homogeneous, else False.
-        inline bool isHomogeneous() const {
-            return mvData.size() < 2; // only one element, or zero element on the list
-        }
-
-        /// \brief inplace simplify the multivector such that all the values with a magnitude lower than a epsilon in the Mv are set to 0.
-        /// \param epsilon - threshold, with default value the epsilon of the float/double/long double type from numeric_limits
-        void roundZero(const T epsilon = std::numeric_limits<T>::epsilon());
-
-        /// \brief Specify if two multivectors have the same grade
-        /// \param mv1 - first multivector
-        /// \param mv2 - second multivector
-        /// \return true if the two multivectors have the same grade, else return false
-        inline bool sameGrade(const Mvec<T>& mv) const{
-            return grade() == mv.grade();
-        }
-
-        /// \brief Display the multivector data (per grade value)
-        void display() const;
-
-        /// \cond DEV
-        /// \brief a function to output all the element of a k-vector in the multivector indexing order.
-        /// \tparam T - type of the multivector
-        /// \param stream - stream that will contain the result
-        /// \param mvec - multivector to be processed
-        /// \param gradeMV - the considered grade
-        /// \param moreThanOne - true if it the first element to display (should we put a '+' before)
-        template<typename U>
-        friend void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1>, unsigned int gradeMV, bool& moreThanOne);
-        /// \endcond // do not comment this functions
-
-/*
-        /// \cond DEV
-        /// \brief functions that enables to  to be able to extract a component of a multivector as a multivector. As an example, mv1.e(1,3) will create a multivector whose 1,3 component will be the component of mv1.
-        /// \tparam Arguments - denotes a variadic list of index
-        /// \param listIndices - the list of indices
-        /// \return a homogeneous multivector filled with zero except the listIndices index
-        /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high
-        template<typename... List>
-        Mvec e(List ...listIndices){
-            // operation on these arguments
-            // First determine the grade, simply use the variadic function sizeof...()
-            const int gd = sizeof...(listIndices); //
-
-            // Second, from these parameters, compute the index in the corresponding VectorXd
-            const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);
-
-            Mvec res;
-            res.mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);
-            res.mvData[gd](idx) = mvData[gd](idx);
-
-            return res;
-        }
-        /// \endcond // do not comment this functions
-*/
-
-
-        /// \brief functions that enables to extract one component of this multivector and initialize another mv with this component
-        /// \tparam Arguments - denotes a variadic list of index
-        /// \param grade : the considered grade
-        /// \param sizeOfKVector: C(dimension, grade)
-        /// \param indexInKvector: the considerd index
-        /// \return a new mv with the right component at the right place
-        Mvec extractOneComponent(const int grade, const int sizeOfKVector, const int indexInKvector) const {
-            Mvec mv;
-            auto itThisMV = this->findGrade(grade);
-            if(itThisMV==this->mvData.end()){
-                return mv;
-            }
-            mv = mv.componentToOne(grade,indexInKvector);
-            mv.mvData.begin()->vec.coeffRef(indexInKvector) = itThisMV->vec.coeff(indexInKvector);
-            return mv;
-        }
-
-
-        /// \brief set of functions that return multivectors containing only the value '1' for the specified component.
-project_multivector_one_component
-
-    };  // end of the class definition
-
-
-    /* ------------------------------------------------------------------------------------------------ */
-
-
-    template<typename T>
-    Mvec<T>::Mvec():gradeBitmap(0)
-    {}
-
-
-    template<typename T>
-    Mvec<T>::Mvec(const Mvec& mv) : mvData(mv.mvData), gradeBitmap(mv.gradeBitmap)
-    {
-        //std::cout << "copy constructor " << std::endl;
-    }
-
-
-    template<typename T>
-    Mvec<T>::Mvec(Mvec<T>&& multivector) : mvData(std::move(multivector.mvData)), gradeBitmap(multivector.gradeBitmap)
-    {
-        // the move constructor
-        //std::cout << "move constructor" << std::endl;
-    }
-
-
-    template<typename T>
-    template<typename U>
-    Mvec<T>::Mvec(const Mvec<U> &mv) : gradeBitmap(mv.gradeBitmap)
-    {
-        for(auto it=mv.mvData.begin(); it != mv.mvData.end(); ++it){
-            Kvec<T> kvec;
-            kvec.grade = it->grade;
-            kvec.vec = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[it->grade]);
-            for(unsigned int i=0; i<it->vec.size(); ++i)
-                kvec.vec.coeffRef(i) = T(it->vec.coeff(i));
-            mvData.push_back(kvec);
-        }
-    }
-
-
-    template<typename T>
-    template<typename U>
-    Mvec<T>::Mvec(const U val) {
-        if(val != U(0)) {
-            gradeBitmap = 1;
-            mvData.push_back({.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>(1),.grade=0});
-            mvData.begin()->vec.coeffRef(0) = val;
-        }
-    }
-
-
-    template<typename T>
-    Mvec<T>::~Mvec()
-    {}
-
-
-    template<typename T>
-    Mvec<T>& Mvec<T>::operator=(const Mvec& mv){
-        if(&mv == this) return *this;
-        gradeBitmap = mv.gradeBitmap;
-        mvData = mv.mvData;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T>& Mvec<T>::operator=(Mvec&& mv){
-        mvData = std::move(mv.mvData);
-        gradeBitmap = mv.gradeBitmap;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator+(const Mvec<T> &mv2) const {
-        Mvec<T> mv3(*this);
-        for(auto & itMv : mv2.mvData) {
-            auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);
-            it->vec += itMv.vec;
-        }
-        return mv3;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator+(const S &value, const Mvec<U> &mv){
-        return mv + value;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator+(const S &value) const {
-        Mvec<T> mv(*this);
-        if(value != T(0)) {
-            auto it = mv.createVectorXdIfDoesNotExist(0);
-            it->vec.coeffRef(0) += value;
-        }
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T>& Mvec<T>::operator+=(const Mvec& mv){
-        *this = *this + mv;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> operator-(const Mvec<T> &mv) { // unary -
-        Mvec<T> mv2(mv);
-        for(auto & itMv : mv2.mvData)
-            itMv.vec = -itMv.vec;
-        return mv2;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator-(const Mvec<T> &mv2) const {
-        Mvec<T> mv3(*this);
-        for(auto & itMv : mv2.mvData) {
-            auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);
-            it->vec -= itMv.vec;
-        }
-        return mv3;
-    }
-
-
-    template<typename T>
-    Mvec<T>& Mvec<T>::operator-=(const Mvec& mv){
-        *this = *this - mv;
-        return *this;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator-(const S &value, const Mvec<U> &mv){
-        return (-mv) + value;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator-(const S &value) const {
-        Mvec<T> mv(*this);
-        if(value != T(0)) {
-            auto it = mv.createVectorXdIfDoesNotExist(0);
-            it->vec.coeffRef(0) -= value;
-        }
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator^(const Mvec<T> &mv2) const {
-#if 0 // only recursive version
-        // Loop over non-empty grade of mv1 and mv2
-        // This version (with recursive call) is only faster than the standard recursive call if
-        // the multivector mv1 and mv2 are homogeneous or near from homogeneous.
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
-            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
-            {
-                unsigned int grade_mv3 = itMv1.grade + itMv2.grade;
-                if(grade_mv3 <= algebraDimension) {
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
-                    outerProductHomogeneous(itMv1.vec, itMv2.vec, itMv3->vec,
-                                            itMv1.grade, itMv2.grade, grade_mv3);
-                }
-            }
-        return mv3;
-#else // use the adaptative pointer function array
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled outer function using the functions pointer called outerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-                if(itMv1.grade + itMv2.grade <= (int) algebraDimension ){
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(itMv1.grade + itMv2.grade);
-                    outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-                }
-            }
-        return mv3;
-#endif
-    }
-
-    template<typename U, typename S>
-    Mvec<U> operator^(const S &value, const Mvec<U> &mv) {
-        return  mv^value;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator^(const S &value) const {
-        Mvec<T> mv2(*this);
-        for(auto & itMv : mv2.mvData)
-            itMv.vec *= T(value);
-        return mv2;
-    }
-
-
-    template<typename T>
-    Mvec<T> &Mvec<T>::operator^=(const Mvec &mv) {
-        *this = *this ^ mv;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator|(const Mvec<T> &mv2) const{
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-
-                // perform the inner product
-                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
-                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
-                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-
-                // check if the result is non-zero
-                if(!((itMv3->vec.array() != 0.0).any())){
-                    mv3.mvData.erase(itMv3);
-                    mv3.gradeBitmap &= ~(1<<absGradeMv3);
-                }
-            }
-        return mv3;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator|(const S &value, const Mvec<U> &mv) {
-        return mv | value;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator|(const S &value) const {
-        return (*this) > value ; // inner between a mv and a scalar gives 0 (empty multivector)
-    }
-
-
-    template<typename T>
-    Mvec<T> &Mvec<T>::operator|=(const Mvec &mv) {
-        *this = *this | mv;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator>(const Mvec<T> &mv2) const{
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-
-                // right contraction constraint: gradeMv1 >= gradeMv2
-                if(itMv1.grade < itMv2.grade)
-                    continue;
-
-                // perform the inner product
-                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
-                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
-                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-
-                // check if the result is non-zero
-                if(!((itMv3->vec.array() != 0.0).any())){
-                    mv3.mvData.erase(itMv3);
-                    mv3.gradeBitmap &= ~(1<<absGradeMv3);
-                }
-            }
-
-        return mv3;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator>(const S &value, const Mvec<U> &mv) {
-        if( (mv.gradeBitmap & 1) == 0) return Mvec<U>();
-        else return Mvec<U>( U(mv.mvData.begin()->vec.coeff(0) * value ) );
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator>(const S &value) const {
-        // return mv x value
-        Mvec<T> mv(*this);
-        for(auto & itMv : mv.mvData)
-            itMv.vec *= T(value);
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator<(const Mvec<T> &mv2) const{
-
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-
-                // left contraction constraint: gradeMv1 <= gradeMv2
-                if(itMv1.grade > itMv2.grade)
-                    continue;
-
-                // perform the inner product
-                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
-                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
-                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-
-                // check if the result is non-zero
-                if(!((itMv3->vec.array() != 0.0).any())){
-                    mv3.mvData.erase(itMv3);
-                    mv3.gradeBitmap &= ~(1<<absGradeMv3);
-                }
-            }
-        return mv3;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator<(const S &value, const Mvec<U> &mv) {
-        // return mv x value
-        Mvec<U> mv3(mv);
-        for(auto & itMv : mv3.mvData)
-            itMv.vec *= U(value);
-        return mv3;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator<(const S &value) const {
-        if( (gradeBitmap & 1) == 0) return Mvec<T>();
-        else return Mvec<T>( T(mvData.begin()->vec.coeff(0) * value ) );
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::scalarProduct(const Mvec<T> &mv2) const{
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-
-                if(itMv1.grade != itMv2.grade)
-                    continue;
-
-                // perform the inner product
-                int absGradeMv3 = 0;
-                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
-                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-
-                // check if the result is non-zero
-                if(!((itMv3->vec.array() != 0.0).any())){
-                    mv3.mvData.erase(itMv3);
-                    mv3.gradeBitmap &= ~(1<<absGradeMv3);
-                }
-            }
-        return mv3;
-    }
-
-
-project_singular_metric_comment_begin
-    template<typename T>
-    Mvec<T> Mvec<T>::outerPrimalDual(const Mvec<T> &mv2) const{
-        Mvec<T> mv3;
-
-        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
-            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
-            {
-                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
-                if(grade_mv3 <= algebraDimension) {
-                    // handle the scalar product as well as the left contraction
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
-                    outerProductPrimalDual(itMv1.vec, itMv2.vec, itMv3->vec,
-                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
-
-                    // check if the result is non-zero
-                    if(!((itMv3->vec.array() != 0.0).any())){
-                        mv3.mvData.erase(itMv3);
-                        mv3.gradeBitmap &= ~(1<<grade_mv3);
-                    }
-                }
-            }
-
-        return mv3;
-    }
-
-    template<typename T>
-    Mvec<T> Mvec<T>::outerDualPrimal(const Mvec<T> &mv2) const{
-        Mvec<T> mv3;
-
-        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
-            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
-            {
-                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
-                if(grade_mv3 <= algebraDimension) {
-                    // handle the scalar product as well as the left contraction
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
-                    outerProductDualPrimal(itMv1.vec, itMv2.vec, itMv3->vec,
-                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
-
-                    // check if the result is non-zero
-                    if(!((itMv3->vec.array() != 0.0).any())){
-                        mv3.mvData.erase(itMv3);
-                        mv3.gradeBitmap &= ~(1<<grade_mv3);
-                    }
-                }
-            }
-
-
-        return mv3;
-
-    }
-
-    template<typename T>
-    Mvec<T> Mvec<T>::outerDualDual(const Mvec<T> &mv2) const{
-        Mvec<T> mv3;
-
-        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1
-            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2
-            {
-                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);
-                if(grade_mv3 <= algebraDimension) {
-                    // handle the scalar product as well as the left contraction
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);
-                    outerProductDualDual(itMv1.vec, itMv2.vec, itMv3->vec,
-                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));
-
-                    // check if the result is non-zero
-                    if(!((itMv3->vec.array() != 0.0).any())){
-                        mv3.mvData.erase(itMv3);
-                        mv3.gradeBitmap &= ~(1<<grade_mv3);
-                    }
-                }
-            }
-
-
-        return mv3;
-
-    }
-
-project_singular_metric_comment_end
-
-    template<typename T>
-    Mvec<T> Mvec<T>::hestenesProduct(const Mvec<T> &mv2) const{
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-
-                if(itMv1.grade*itMv2.grade == 0)
-                    continue;
-
-                // perform the inner product
-                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));
-                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);
-                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-
-                // check if the result is non-zero
-                if(!((itMv3->vec.array() != 0.0).any())){
-                    mv3.mvData.erase(itMv3);
-                    mv3.gradeBitmap &= ~(1<<absGradeMv3);
-                }
-            }
-        return mv3;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator*(const Mvec<T> &mv2) const {
-        // Loop over non-empty grade of mv1 and mv2
-        // call the right explicit unrolled product function using the functions pointer
-        Mvec<T> mv3;
-        for(const auto & itMv1 : this->mvData)
-            for(const auto & itMv2 : mv2.mvData){
-                project_select_recursive_geometric_product_template
-                // outer product block
-                unsigned int gradeOuter = itMv1.grade + itMv2.grade;
-                if(gradeOuter <=  algebraDimension ){
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeOuter);
-                    outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-                    if(!((itMv3->vec.array() != 0.0).any())){
-                        mv3.mvData.erase(itMv3);
-                        mv3.gradeBitmap &= ~(1<<gradeOuter);
-                    }
-                }
-
-                // inner product block
-                unsigned int gradeInner = (unsigned int)std::abs((int)(itMv1.grade-itMv2.grade));
-                // when the grade of one of the kvectors is zero, the inner product is the same as the outer product
-                if(gradeInner != gradeOuter) {
-                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeInner);
-                    innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-                    // check if the result is non-zero
-                    if(!((itMv3->vec.array() != 0.0).any())){
-                        mv3.mvData.erase(itMv3);
-                        mv3.gradeBitmap &= ~(1<<gradeInner);
-                    }
-
-                    // geometric product part
-                    int gradeMax = std::min(((2*algebraDimension)-gradeOuter)+1,gradeOuter);
-                    for (int gradeResult = gradeInner+2; gradeResult < gradeMax; gradeResult+=2) {
-                        auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeResult);
-                        geometricFunctionsContainer<T>[gradeResult][itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);
-                        // check if the result is non-zero
-                        if(!((itMv3->vec.array() != 0.0).any())){
-                            mv3.mvData.erase(itMv3);
-                            mv3.gradeBitmap &= ~(1<<gradeResult);
-                        }
-                    }
-                }
-            }
-        return mv3;
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator*(const S &value, const Mvec<U> &mv){
-        return mv * value;
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator*(const S &value) const{
-        Mvec<T> mv(*this);
-        for(auto & itMv : mv.mvData)
-            itMv.vec *= T(value);
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T> &Mvec<T>::operator*=(const Mvec &mv) {
-        *this = *this * mv;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::operator/(const Mvec<T> &mv2) const {
-        return *this * mv2.inv();
-    }
-
-
-    template<typename U, typename S>
-    Mvec<U> operator/(const S &value, const Mvec<U> &mv){
-        return value * mv.inv();
-    }
-
-
-    template<typename T>
-    template<typename S>
-    Mvec<T> Mvec<T>::operator/(const S &value) const {
-        Mvec<T> mv(*this);
-        for(auto & itMv : mv.mvData)
-            itMv.vec /= value;
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T> &Mvec<T>::operator/=(const Mvec &mv) {
-        *this = *this / mv;
-        return *this;
-    }
-
-
-    template<typename T>
-    Mvec<T> operator~(const Mvec<T> &mv){
-        return mv.reverse();
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::inv() const {
-        T n = this->quadraticNorm();
-        if(n<std::numeric_limits<T>::epsilon() && n>-std::numeric_limits<T>::epsilon())
-            return Mvec<T>(); // return 0, this is was gaviewer does.
-        return this->reverse() / n;
-    };
-
-
-    project_singular_metric_comment_begin
-    template<typename U>
-    Mvec<U> operator!(const Mvec<U> &mv) {
-        return mv.dual();
-    }
-project_singular_metric_comment_end
-
-    /// \brief defines the left contraction between two multivectors
-    /// \param mv1 - a multivector
-    /// \param mv2 - second operand, corresponds to a multivector
-    /// \return mv1 left contraction mv2
-    template<typename T>
-    Mvec<T> leftContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){
-        return mv1 < mv2;
-    }
-
-
-    /// \brief defines the right contraction between two multivectors
-    /// \param mv1 - a multivector
-    /// \param mv2 - second operand, corresponds to a multivector
-    /// \return mv1 right contraction mv2
-    template<typename T>
-    Mvec<T> rightContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){
-        return mv1 > mv2;
-    }
-
-
-    /// \brief returns a multivector that only contains the coefficient associated to the pseudoscalar.
-    /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.
-    template<typename T>
-    Mvec<T> I() {
-        Mvec<T> mvec;
-        mvec[project_pseudo_scalar] = 1.0;
-        return mvec;
-    }
-
-project_singular_metric_comment_begin
-    /// \brief return the inverse of the pseudo scalar.
-    /// \return returns a multivector corresponding to the inverse of the pseudo scalar.
-    template<typename T>
-    Mvec<T> Iinv() {
-        Mvec<T> mvec;
-        mvec[project_pseudo_scalar] = pseudoScalarInverse; // we do not return just a scalar of type T because the grade should be dimension and not 0.
-        return mvec;
-    }
-project_singular_metric_comment_end
-
-project_singular_metric_comment_begin
-    // compute the dual of a multivector (i.e mv* = mv.reverse() * Iinv)
-    template<typename T>
-    Mvec<T> Mvec<T>::dual() const {
-        Mvec<T> mvResult;
-        // for each k-vectors of the multivector
-        for(auto itMv=mvData.rbegin(); itMv!=mvData.rend(); ++itMv){
-
-            // create the dual k-vector
-            Kvec<T> kvec ={itMv->vec,algebraDimension-(itMv->grade)};
-
-            // some elements need to be permuted
-            for(unsigned int i=0;i<binomialArray[itMv->grade];++i)
-                kvec.vec.coeffRef(dualPermutations[itMv->grade][i]) = itMv->vec.coeff(i);
-
-            // the inner product may involve some constant multiplucation for the dual elements
-            kvec.vec = kvec.vec.cwiseProduct(dualCoefficients[itMv->grade].template cast<T>());
-
-            // add the k-vector to the resulting multivector
-            mvResult.mvData.push_back(kvec);
-            mvResult.gradeBitmap |= (1 << kvec.grade);
-        }
-        return mvResult;
-    };
-project_singular_metric_comment_end
-
-    // \brief compute the reverse of a multivector
-    // \return - the reverse of the multivector
-    template<typename T>
-    Mvec<T> Mvec<T>::reverse() const {
-        Mvec<T> mv(*this);
-        for(auto & itMv : mv.mvData)
-            if(signReversePerGrade[itMv.grade] == -1)
-                itMv.vec *= -1;
-        return mv;
-    }
-
-
-    template<typename T>
-    Mvec<T> Mvec<T>::grade(const int i) const{
-
-        auto it = findGrade(i);
-        Mvec<T> mv;
-
-        // if not found, return the empty multivector
-        if(it == mvData.end())
-            return mv;
-
-        // else return the grade 'i' data
-        mv.mvData.push_back(*it);
-        mv.gradeBitmap = 1 << (i);
-
-        return mv;
-    }
-
-
-    template<typename T>
-    std::vector<unsigned int> Mvec<T>::grades() const {
-        std::vector<unsigned int> mvGrades;
-        for(unsigned int i=0; i< algebraDimension+1; ++i){
-
-            if( (gradeBitmap & (1<<i)))
-                mvGrades.push_back(i);
-        }
-
-        return mvGrades;
-    }
-
-
-    template<typename T>
-    void Mvec<T>::roundZero(const T epsilon) {
-        // loop over each k-vector of the multivector
-        auto itMv=mvData.begin();
-        while(itMv != mvData.end()){
-            // loop over each element of the k-vector
-            for(unsigned int i=0; i<(unsigned int)itMv->vec.size(); ++i)
-                if(fabs(itMv->vec.coeff(i)) <= epsilon)
-                    itMv->vec.coeffRef(i) = 0.0;
-            // if the k-vector is full of 0, remove it
-            if(!((itMv->vec.array() != 0.0).any())){
-                gradeBitmap = gradeBitmap - (1 << itMv->grade);
-                mvData.erase(itMv++);
-            }
-            else ++itMv;
-        }
-    }
-
-
-    template<typename T>
-    void Mvec<T>::clear(const int grade) {
-
-        // full erase
-        if(grade < 0){
-            mvData.clear();
-            gradeBitmap = 0;
-            return;
-        }
-
-        // partial erase
-        auto iter = findGrade(grade);
-        if(iter != mvData.end()) {
-            gradeBitmap = gradeBitmap - (1 << iter->grade);
-            iter = mvData.erase(iter);
-        }
-    }
-
-
-    template<typename T>
-    void Mvec<T>::display() const {
-        if(mvData.size() == 0)
-            std::cout << " null grade , null value " <<std::endl;
-
-        for(auto itmv=mvData.begin(); itmv!=mvData.end(); ++itmv){
-            std::cout << "  grade   : " << itmv->grade  << std::endl;
-            std::cout << "  kvector : " << itmv->vec.transpose() << std::endl;
-        }
-        std::cout << std::endl;
-    }
-
-
-    /// \cond DEV
-    template<typename U>
-    void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1> kvector, unsigned int gradeMV, bool& moreThanOne ){
-
-        // version with XOR indices
-        // for the current grade, generate all the XOR indices
-        std::vector<bool> booleanXorIndex(algebraDimension);
-        std::fill(booleanXorIndex.begin(), booleanXorIndex.end(), false);
-        // build the first xorIndex of grade 'grade' (with 'grade' values to true).
-        std::fill(booleanXorIndex.begin(), booleanXorIndex.begin() + gradeMV, true);
-        unsigned int positionInKVector = 0;
-
-        do {
-            // convert the vector of bool to the string containing all the basis vectors
-            std::string basisBlades={};
-            for(unsigned int i=0; i<algebraDimension; ++i) {
-                if(booleanXorIndex[i]) {
-                    basisBlades += basisVectors[i];
-                }
-            }
-
-            if(kvector.coeff(positionInKVector)!= 0) {
-                if(!(moreThanOne)){
-                    stream<< kvector.coeff(positionInKVector) << "*e"+ basisBlades;
-                    moreThanOne = true;
-                }else{
-                    if(kvector.coeff(positionInKVector)>0)
-                        stream<< " + " << kvector.coeff(positionInKVector) << "*e" + basisBlades;
-                    else
-                        stream<< " - " << -kvector.coeff(positionInKVector) << "*e" + basisBlades;
-                }
-            }
-
-            positionInKVector++;
-
-        } while(std::prev_permutation(booleanXorIndex.begin(), booleanXorIndex.end())); // compute next permutation of the true values of booleanXorIndex
-
-    }
-    /// \endcond
-
-
-    template<typename U>
-    std::ostream &operator<<(std::ostream &stream, const Mvec<U> &mvec) {
-        if(mvec.mvData.size()==0){
-            stream << " 0 ";
-            return stream;
-        }
-
-        bool moreThanOne = false;
-        auto mvIterator = mvec.mvData.begin();
-
-        // if the iterator corresponds to the scalar
-        if(mvIterator->grade == 0){
-            stream << mvIterator->vec.coeff(0);
-            moreThanOne = true;
-            ++mvIterator;
-        }
-
-        // for all other k-vectors of mvec
-        for(; mvIterator!=mvec.mvData.end(); ++mvIterator){
-
-            // call the function that covers a single k-vector
-            traverseKVector(stream,mvIterator->vec,mvIterator->grade, moreThanOne);
-        }
-
-        if(!moreThanOne)
-            stream << " 0 ";
-
-        return stream;
-    }
-
-    // static unit multivector functions (returns a multivector with only one component)
-project_static_multivector_one_component
-
-    //    template<typename U>
-    //    void recursiveTraversalMultivector(std::ostream &stream, const Mvec<U> &mvec, unsigned int currentGrade, int currentIndex, std::vector<int> listBasisBlades, unsigned int lastIndex, unsigned int gradeMV, bool& moreThanOne);
-
-
-    void temporaryFunction1();
-
-}     /// End of Namespace
-
-#endif // project_inclusion_guard
+// Copyright (c) 2018 by University Paris-Est Marne-la-Vallee\r
+// Mvec.hpp\r
+// This file is part of the Garamon for project_namespace.\r
+// Authors: Stephane Breuils and Vincent Nozick\r
+// Contact: vincent.nozick@u-pem.fr\r
+//\r
+// Licence MIT\r
+// A a copy of the MIT License is given along with this program\r
+\r
+/// \file Mvec.hpp\r
+/// \author Stephane Breuils, Vincent Nozick\r
+/// \brief Class to define a multivector and its basic operators in the Geometric algebra of project_namespace.\r
+\r
+\r
+// Anti-doublon\r
+#ifndef project_inclusion_guard\r
+#define project_inclusion_guard\r
+#pragma once\r
+\r
+// External Includes\r
+#include <Eigen/Core>\r
+#include <list>\r
+#include <iostream>\r
+#include <cmath>\r
+#include <limits>\r
+\r
+// Internal Includes\r
+#include "project_namespace/Utility.hpp"\r
+#include "project_namespace/Constants.hpp"\r
+\r
+#include "project_namespace/Outer.hpp"\r
+#include "project_namespace/Inner.hpp"\r
+#include "project_namespace/Geometric.hpp"\r
+\r
+#include "project_namespace/OuterExplicit.hpp"\r
+#include "project_namespace/InnerExplicit.hpp"\r
+#include "project_namespace/GeometricExplicit.hpp"\r
+\r
+/*!\r
+ * @namespace project_namespace\r
+ */\r
+namespace project_namespace{\r
+\r
+\r
+    /// \class Kvec\r
+    /// \brief class defining a single grade component of a multivector.\r
+    template<class T>\r
+    struct Kvec{\r
+\r
+        Eigen::Matrix<T, Eigen::Dynamic, 1> vec;  /*!< dynamic vector of Eigen Library */\r
+\r
+        unsigned int grade; /*!< grade k of the k-vector */\r
+\r
+        /// \brief operator == to test the equality between 2 k-vectors.\r
+        bool operator == (const Kvec<T>& other) const {\r
+            return vec == other.vec;\r
+        }\r
+    };\r
+\r
+\r
+\r
+    /// \class Mvec\r
+    /// \brief class defining multivectors.\r
+    template<typename T = double>\r
+    class Mvec {\r
+\r
+    protected:\r
+        std::list<Kvec<T>> mvData;  /*!< set of k-vectors, mapped by grade */\r
+        unsigned int gradeBitmap;   /*!< ith bit to 1 if grade i is contained in the multivector */\r
+\r
+    public:\r
+\r
+        /// \brief Default constructor, generate an empty multivector equivalent to the scalar 0.\r
+        Mvec();\r
+\r
+        /// \brief Copy constructor\r
+        /// \param mv - the multivector to be copied\r
+        Mvec(const Mvec& mv);\r
+\r
+        /// \brief Copy constructor of Mvec\r
+        /// \param mv - the multivector which has to be copied\r
+        Mvec(Mvec&& mv); // move constructor\r
+\r
+\r
+        template <typename U>\r
+        friend class Mvec;\r
+\r
+        /// \brief Copy constructor of Mvec from different types\r
+        /// \param mv - the multivector with another template type\r
+        template<typename U>\r
+        Mvec<T>(const Mvec<U> &mv);\r
+\r
+        /// \brief Constructor of Mvec from a scalar\r
+        /// \param val - scalar value\r
+        template<typename S>\r
+        Mvec(const S val);\r
+\r
+        /// Destructor\r
+        ~Mvec();\r
+\r
+        /// \brief Overload the assignment operator. No need any copy when the argument is an R-value, just need to move.\r
+        /// \param mv - Mvec used as a Rvalue\r
+        /// \return assign other to this\r
+        Mvec& operator=(Mvec&& mv);\r
+\r
+        /// \brief Overload the assignment operator\r
+        /// \param mv - Mvec\r
+        /// \return assign mv to this object\r
+        Mvec& operator=(const Mvec& mv);\r
+\r
+        /// \brief defines the addition between two Mvec\r
+        /// \param mv2 - second operand of type Mvec\r
+        /// \return this + mv2\r
+        Mvec operator+(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the addition between a Mvec and a scalar\r
+        /// \param value - second operand (scalar)\r
+        /// \return this + scalar\r
+        template<typename S>\r
+        Mvec operator+(const S &value) const;\r
+\r
+        /// \brief defines the addition between a scalar and a Mvec\r
+        /// \param value a scalar\r
+        /// \param mv - second operand of type Mvec\r
+        /// \return scalar + mv\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator+(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the += operator, corresponds to this += mv\r
+        /// \param mv - Mvec to be added to this object\r
+        /// \return this += mv\r
+        Mvec& operator+=(const Mvec& mv);\r
+\r
+        /// \brief defines the opposite of a multivector\r
+        /// \param mv: operand of type Mvec\r
+        /// \return -mv\r
+        template<typename U>\r
+        friend Mvec<U> operator-(const Mvec<U> &mv); // unary operator -mv1\r
+\r
+        /// \brief defines the difference between two Mvec\r
+        /// \param mv2 - second operand of type Mvec\r
+        /// \return this - mv2\r
+        Mvec<T> operator-(const Mvec<T> &mv2) const;\r
+\r
+        /// \brief defines the difference between a Mvec and a scalar\r
+        /// \param value - second operand (scalar)\r
+        /// \return this - scalar\r
+        template<typename S>\r
+        Mvec operator-(const S &value) const;\r
+\r
+        /// \brief defines the difference between a scalar and a Mvec\r
+        /// \param value a scalar\r
+        /// \param mv - second operand of type Mvec\r
+        /// \return scalar - mvX\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator-(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the -= operator, corresponds to this -= mv\r
+        /// \param mv - Mvec to be added to this object\r
+        /// \return this -= mv\r
+        Mvec& operator-=(const Mvec& mv);\r
+\r
+        /// \brief defines the outer product between two multivectors\r
+        /// \param mv2 - a multivector\r
+        /// \return this^mv2\r
+        Mvec operator^(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the outer product a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return this^value\r
+        template<typename S>\r
+        Mvec operator^(const S &value) const;\r
+\r
+        /// \brief defines the outer product between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return value ^ mv\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator^(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the outer product with operator =, corresponds to this ^= mv\r
+        /// \param mv - Mvec to be wedged to this object\r
+        /// \return this ^= mv\r
+        Mvec& operator^=(const Mvec& mv);\r
+\r
+        /// \brief defines the inner product between two multivectors\r
+        /// \param mv2 - a multivector\r
+        /// \return this.mv2\r
+        Mvec operator|(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the inner product between a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return this.value = 0 (empty multivector)\r
+        template<typename S>\r
+        Mvec operator|(const S &value) const;\r
+\r
+        /// \brief defines the inner product between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return mv.value = 0 (empty multivector)\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator|(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the inner product with operator =, corresponds to this |= mv\r
+        /// \param mv - Mvec to be computed in the inner product with this object\r
+        /// \return this |= mv\r
+        Mvec& operator|=(const Mvec& mv);\r
+\r
+        /// \brief defines the right contraction between two multivectors\r
+        /// \param mv2 - a multivector\r
+        /// \return the right contraction : $this \\lfloor mv2$\r
+        Mvec operator>(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the right contraction between a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return $this \\lfloor value$\r
+        template<typename S>\r
+        Mvec operator>(const S &value) const;\r
+\r
+        /// \brief defines the right contraction between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return $value \\lfloor this = 0$ (empty multivector)\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator>(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief defines the left contraction between two multivectors\r
+        /// \param mv2 - a multivector\r
+        /// \return the left contraction $this \\rfloor mv2$\r
+        Mvec operator<(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the left contraction between a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return $value \\rfloor this$ = 0 (empty multivector)\r
+        template<typename S>\r
+        Mvec operator<(const S &value) const;\r
+\r
+        /// \brief defines the left contraction between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return $this \\rfloor value$\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator<(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief defines the geometric product between two multivectors\r
+        /// \param mv2 - a multivector\r
+        /// \return mv1*mv2\r
+        Mvec operator*(const Mvec &mv2) const;\r
+\r
+    project_singular_metric_comment_begin\r
+        /// \brief defines the outer product between two multivectors, where the second multivector is dualized during the product computation.\r
+        /// \param mv2 - a multivector that will be dualized during the wedge with the calling multivector.\r
+        /// \return a multivector.\r
+        Mvec<T> outerPrimalDual(const Mvec<T> &mv2) const;\r
+\r
+        /// \brief defines the outer product between two multivectors, where the first multivector is dualized during the product computation.\r
+        /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.\r
+        /// \return a multivector.\r
+        Mvec<T> outerDualPrimal(const Mvec<T> &mv2) const;\r
+\r
+        /// \brief defines the outer product between two multivectors, where both multivectors are dualized during the product computation.\r
+        /// \param mv2 - a primal form of a multivector; the object multivector will be dualized during the wedge with the calling multivector.\r
+        /// \return a multivector.\r
+        Mvec<T> outerDualDual(const Mvec<T> &mv2) const;\r
+\r
+    project_singular_metric_comment_end\r
+        \r
+\r
+        /// \brief defines the scalar product between two multivectors (sum of the inner products between same grade pairs from the 2 multivectors)\r
+        /// \param mv2 - a multivector\r
+        /// \return a scalar\r
+        Mvec<T> scalarProduct(const Mvec<T> &mv2) const;\r
+\r
+        /// \brief defines the Hestenes product between two multivectors (Inner product - scalar product)\r
+        /// \param mv2 - a multivector\r
+        /// \return a multivector.\r
+        Mvec<T> hestenesProduct(const Mvec<T> &mv2) const;\r
+\r
+        /// \brief defines the geometric product between a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return mv2*value\r
+        template<typename S>\r
+        Mvec operator*(const S &value) const;\r
+\r
+        /// \brief defines the geometric product between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return value*mv2\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator*(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the geometric product with operator =, corresponds to this *= mv\r
+        /// \param mv - Mvec to be multiplied to this object\r
+        /// \return this *= mv\r
+        Mvec& operator*=(const Mvec& mv);\r
+\r
+        /// \brief defines the geometric product with a multivector and the inverse of a second multivector\r
+        /// \param mv2 - a multivector\r
+        /// \return this / mv2\r
+        Mvec operator/(const Mvec &mv2) const;\r
+\r
+        /// \brief defines the scalar inverse between a multivector and a scalar\r
+        /// \param value - a scalar\r
+        /// \return this / value\r
+        template<typename S>\r
+        Mvec operator/(const S &value) const;\r
+\r
+        /// \brief defines the inverse product between a scalar and a multivector\r
+        /// \param value - a scalar\r
+        /// \param mv - a multivector\r
+        /// \return value / mv\r
+        template<typename U, typename S>\r
+        friend Mvec<U> operator/(const S &value, const Mvec<U> &mv);\r
+\r
+        /// \brief Overload the inverse with operator =, corresponds to this /= mv\r
+        /// \param mv - Mvec to be inversed to this object\r
+        /// \return this /= mv\r
+        Mvec& operator/=(const Mvec& mv);\r
+\r
+        /// \brief the reverse of a multivector, i.e. if mv = a1^a2^...^an, then reverse(mv) = an^...^a2^a1\r
+        /// \param mv - a multivector\r
+        /// \return reverse of mv\r
+        template<typename U>\r
+        friend Mvec<U> operator~(const Mvec<U> &mv);\r
+\r
+    project_singular_metric_comment_begin\r
+        /// \brief the dual of a k-vector is defined as $A_k^* = A_k \\lcont I_n^{-1}$, for a multivector, we just dualize all its components.\r
+        /// \param mv - a multivector\r
+        /// \return the dual of mv\r
+        template<typename U>\r
+        friend Mvec<U> operator!(const Mvec<U> &mv);\r
+    project_singular_metric_comment_end\r
+        /// \brief boolean operator that tests the equality between two Mvec\r
+        /// \param mv2 - second operand of type Mvec\r
+        /// \return whether two Mvec have the same coefficients\r
+        inline bool operator==(const Mvec& mv2){\r
+            if(gradeBitmap != mv2.gradeBitmap)\r
+                return false;\r
+\r
+            return mvData == mv2.mvData;   //// listcaca : ca marche que si les listes sont ordonnees\r
+        }\r
+\r
+        /// \brief compute the inverse of a multivector\r
+        /// \return - the inverse of the current multivector\r
+        Mvec<T> inv() const;\r
+\r
+\r
+            /// \brief operator to test whether two Mvec have not the same coefficients\r
+        /// \param mv2 - second operand of type Mvec\r
+        /// \return boolean that specify the non-equality between two Mvec\r
+        inline bool operator!=(const Mvec& mv2){\r
+            return !(mvData == mv2.mvData);\r
+        }\r
+\r
+        /// \brief Display all the non-null basis blades of this objects\r
+        /// \param stream - destination stream\r
+        /// \param mvec - the multivector to be outputed\r
+        /// \return a stream that contains the list of the non-zero element of the multivector\r
+        template<typename U>\r
+        friend std::ostream& operator<< (std::ostream& stream, const Mvec<U> &mvec);\r
+\r
+        /// \brief overload the casting operator, using this function, it is now possible to compute : float a = float(mv);\r
+        /// \return the scalar part of the multivector\r
+        operator T () {\r
+            // if no scalar part, return\r
+\r
+            if( (gradeBitmap & 1) == 0 )\r
+                return 0;\r
+\r
+            // assuming now that the scalar part exists, return it\r
+            return findGrade(0)->vec[0];\r
+        }\r
+\r
+        /// \brief Overload the [] operator to assign a basis blade to a multivector. As an example, float a = mv[E12] = 42.\r
+        /// \param idx - the basis vector index related to the query.\r
+        /// \return the coefficient of the multivector corresponding to the "idx" component.\r
+        /// \todo handle the cases when the number of indices is higher than the dimension, and when the index is too high\r
+        T& operator[](const int idx ){\r
+\r
+            const unsigned int grade = xorIndexToGrade[idx];\r
+            const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];\r
+\r
+            auto it = mvData.begin();\r
+            while(it != mvData.end()){\r
+\r
+                // if grade not reach yet, continue\r
+                if(it->grade < grade)\r
+                    ++it;\r
+                else{\r
+\r
+                    // if grade found\r
+                    if(it->grade == grade)\r
+                        return it->vec[idxHomogeneous];\r
+\r
+                    // if grade exceed, create it and inster it before the current element\r
+                    Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};\r
+                    auto it2 = mvData.insert(it,kvec);\r
+                    gradeBitmap |= 1 << (grade);\r
+                    return it2->vec[idxHomogeneous];\r
+                }\r
+            }\r
+\r
+            // if the searched element should be added at the end, add it\r
+            Kvec<T> kvec = {Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]),grade};\r
+\r
+            auto it2 = mvData.insert(mvData.end(),kvec);\r
+            gradeBitmap |= 1 << (grade);\r
+            return it2->vec[idxHomogeneous];\r
+        }\r
+\r
+        /// \brief Overload the [] operator to copy a basis blade of this multivector. As an example, float a = mv[E12].\r
+        /// \param idx - the basis vector index related to the query.\r
+        /// \return the coefficient of the multivector corresponding to the "idx" component.\r
+        const T& operator[](const int idx) const{\r
+\r
+            const unsigned int grade = xorIndexToGrade[idx];\r
+            const unsigned int idxHomogeneous = xorIndexToHomogeneousIndex[idx];\r
+\r
+            auto it = mvData.begin();\r
+            while(it != mvData.end()){\r
+\r
+                // if grade not reach yet, continue\r
+                if(it->grade < grade)\r
+                    ++it;\r
+                else{\r
+\r
+                    // if grade found\r
+                    if(it->grade == grade)\r
+                        return it->vec[idxHomogeneous];\r
+\r
+                    // if grade exceed, return a reference on ... humm ...\r
+                    T val = T(0);\r
+                    return val;\r
+                }\r
+            }\r
+\r
+            // if the searched element should be added at the end, humm ...\r
+            T val = T(0);   /// \todo find an elegant way to correct this\r
+            return val;\r
+        }\r
+\r
+/*\r
+        /// \cond DEV\r
+        // Variadic operators\r
+        /// \brief Overload the () operator to assign a basis blade to a multivector. As an example, consider we have a Mvec a and we want to set its e23 component to 4.7, using operator(), it will consists in writing a[2,3]=4.8.\r
+        /// \tparam Arguments - denotes a variadic list of index\r
+        /// \param listIndices - the list of indices\r
+        /// \return the right index in the homogeneous vectorXd\r
+        /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high\r
+        template<typename... List>\r
+        T& operator()(List ...listIndices){\r
+            // operation on these arguments\r
+            // First determine the grade, simply use the variadic function sizeof...()\r
+            // This function is correct in terms of indexing\r
+            const int gd = sizeof...(listIndices); //\r
+\r
+            // Second, from these parameters, compute the index in the corresponding VectorXd\r
+            const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);//Binomial<algebraDimension,binomCoef>::value;//binomCoef - sumB(first,listIndices...);//computeOrderingIndex<>(listIndices...);//(Binomial<algebraDimension,gd>::value-1);//-computeOrderingIndex<int,List>(4,2,listIndices...);//computeOrderingIndex<algebraDimension,gd>(listIndices...);//idxVariadic<algebraDimension,gd,first,listIndices...>::value; //computeOrderingIndex<algebraDimension,gd,listIndices...>();//idxVariadic<algebraDimension,gd,listIndices...>::value; // (Binomial<algebraDimension,gd>::value-1)- idxVariadic<algebraDimension,gd,listIndices...>::value\r
+\r
+            if(mvData.count(gd)>0)\r
+                return (mvData.at(gd))(idx);\r
+\r
+            mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);\r
+\r
+            return mvData[gd](idx);\r
+        }\r
+        /// \endcond\r
+*/\r
+\r
+        /// \cond DEV\r
+        /// \brief returns a multivector with one component to 1\r
+        /// \param grade : grade of the component to enable\r
+        /// \param index : index of the parameter in the k-vector (k = grade)\r
+        inline Mvec componentToOne(const unsigned int grade, const int index){\r
+            Kvec<T> kvec;\r
+                       kvec.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]);\r
+                       kvec.grade=grade;\r
+            kvec.vec[index] = T(1);\r
+            Mvec mv1;\r
+            mv1.mvData.push_back(kvec);\r
+            mv1.gradeBitmap = 1 << (grade);\r
+\r
+            return mv1;\r
+        }\r
+        /// \endcond // do not comment this functions\r
+\r
+        /// \cond DEV\r
+        /// \brief create a VectorXd if it has not yet been created\r
+        /// \param grade - grade of the considered kvector\r
+        /// \return nothing\r
+        inline typename std::list<Kvec<T>>::iterator createVectorXdIfDoesNotExist(const unsigned int grade){\r
+\r
+            auto it = mvData.begin();\r
+            while(it != mvData.end()){\r
+\r
+                // if grade not reach yet, continue\r
+                if(it->grade < grade)\r
+                    ++it;\r
+                else{\r
+\r
+                    // if grade found\r
+                    if(it->grade == grade)\r
+                        return it;\r
+\r
+                    // if grade exceed, create it and inster it before the current element\r
+                    Kvec<T> kvec;\r
+                                       kvec.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]);\r
+                                       kvec.grade=grade;\r
+                    auto it2 = mvData.insert(it,kvec);\r
+                    gradeBitmap |= 1 << (grade);\r
+                    return it2;\r
+                }\r
+            }\r
+\r
+            // if the searched element should be added at the end, add it\r
+            Kvec<T> kvec;\r
+                       kvec.vec=Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[grade]);\r
+                       kvec.grade=grade;\r
+            auto it2 = mvData.insert(mvData.end(),kvec);\r
+            gradeBitmap |= 1 << (grade);\r
+            return it2;\r
+        }\r
+\r
+\r
+        /// \cond DEV\r
+        /// \brief modify the element of the multivector whose grade is "grade"\r
+        /// \param grade : the grade to access\r
+        /// \param indexVectorXd : index of the k-vector (with k = "grade")\r
+        /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd\r
+        inline T& at(const int grade, const int indexVectorXd){\r
+            auto it = createVectorXdIfDoesNotExist(grade);\r
+            return it->vec.coeffRef(indexVectorXd);  \r
+        }\r
+        /// \endcond // do not comment this functions\r
+\r
+        /// \cond DEV\r
+        /// \brief return the element of the multivector whose grade is "grade"\r
+        /// \param grade : the grade to access\r
+        /// \param indexVectorXd : index of the k-vector (with k = "grade")\r
+        /// \return the element of the Mv whose grade is grade and index in the VectorXd is indexVectorXd\r
+        inline T at(const int grade, const int indexVectorXd) const{\r
+\r
+            if((gradeBitmap & (1<<(grade))) == 0 )\r
+                return T(0);\r
+\r
+            auto it = findGrade(grade);\r
+            return it->vec.coeff(indexVectorXd);\r
+        }\r
+        /// \endcond // do not comment this functions\r
+\r
+        /// \brief the L2-norm of the mv is sqrt( abs( mv.mv ) )\r
+        /// \return the L2-norm of the multivector (as a double)\r
+        T inline norm() const {\r
+            return sqrt( fabs( (*this).scalarProduct( this->reverse() ) ));\r
+        };\r
+\r
+        /// \brief the L2-norm over 2 of the mv is mv.mv\r
+        /// \return the L2-norm of the multivector (as a double)\r
+        T inline quadraticNorm() const {\r
+            return ( this->reverse() | (*this) );\r
+        };\r
+\r
+    project_singular_metric_comment_begin\r
+        /// \brief compute the dual of a multivector (i.e mv* = reverse(mv) * Iinv)\r
+        /// \return - the dual of the multivector\r
+        Mvec<T> dual() const;\r
+    project_singular_metric_comment_end\r
+        /// \brief compute the reverse of a multivector\r
+        /// \return - the reverse of the multivector\r
+        Mvec<T> reverse() const;\r
+\r
+        /// \brief search in the multivector for a Kvec of grade "grade"\r
+        /// \return return a const iterator on the Kvec if exist, else return mvData.end()\r
+        inline typename std::list<Kvec<T>>::const_iterator findGrade(const unsigned int & gradeToFind) const {\r
+            for(auto it = mvData.begin();it != mvData.end();++it){\r
+                if(it->grade==gradeToFind){\r
+                    return it;\r
+                }\r
+            }\r
+            return mvData.end();\r
+        }\r
+\r
+        /// \brief search in the multivector for a Kvec of grade "grade"\r
+        /// \return return an iterator on the Kvec if exist, else return mvData.end()\r
+        inline typename std::list<Kvec<T>>::iterator findGrade(const unsigned int & gradeToFind) {\r
+            for(auto it = mvData.begin();it != mvData.end();++it){\r
+                if(it->grade==gradeToFind){\r
+                    return it;\r
+                }\r
+            }\r
+            return mvData.end();\r
+        }\r
+\r
+        /// \brief return the (highest) grade of the multivector\r
+        /// \return the highest grade of the multivector\r
+        inline int grade() const {\r
+            return (mvData.begin() == mvData.end()) ? 0 : mvData.rbegin()->grade;\r
+        }\r
+\r
+        /// \brief return the all non-zero grades of the multivector (several grades for non-homogeneous multivectors)\r
+        /// \return a list of the grades encountered in the multivector\r
+        std::vector<unsigned int> grades() const;\r
+\r
+        /// \brief returns a multivector that contains all the components of this multivector whose grade is i\r
+        /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.\r
+        Mvec grade(const int i) const;\r
+\r
+        /// \brief tell whether the multivector has grade component\r
+        /// \param grade - grade of the considered kvector\r
+        /// \return true if multivector has grade component, false else\r
+        inline bool isGrade(const unsigned int grade) const{\r
+            return ( (gradeBitmap & (1<<(grade+1))) != 0 );\r
+        }\r
+\r
+        /// \brief partially or completely erase the content of a multivector\r
+        /// \param grade : if < 0, erase all the multivector, else just erase the part of grade "grade".\r
+        void clear(const int grade = -1);\r
+\r
+        /// \brief check is a mutivector is empty, i.e. corresponds to 0.\r
+        /// \return True if the multivector is empty, else False.\r
+        inline bool isEmpty() const {\r
+            return mvData.empty();\r
+        }\r
+\r
+        /// \brief A multivector is homogeneous if all its components have the same grade.\r
+        /// \return True if the multivector is homogeneous, else False.\r
+        inline bool isHomogeneous() const {\r
+            return mvData.size() < 2; // only one element, or zero element on the list\r
+        }\r
+\r
+        /// \brief inplace simplify the multivector such that all the values with a magnitude lower than a epsilon in the Mv are set to 0.\r
+        /// \param epsilon - threshold, with default value the epsilon of the float/double/long double type from numeric_limits\r
+        void roundZero(const T epsilon = std::numeric_limits<T>::epsilon());\r
+\r
+        /// \brief Specify if two multivectors have the same grade\r
+        /// \param mv1 - first multivector\r
+        /// \param mv2 - second multivector\r
+        /// \return true if the two multivectors have the same grade, else return false\r
+        inline bool sameGrade(const Mvec<T>& mv) const{\r
+            return grade() == mv.grade();\r
+        }\r
+\r
+        /// \brief Display the multivector data (per grade value)\r
+        void display() const;\r
+\r
+        /// \cond DEV\r
+        /// \brief a function to output all the element of a k-vector in the multivector indexing order.\r
+        /// \tparam T - type of the multivector\r
+        /// \param stream - stream that will contain the result\r
+        /// \param mvec - multivector to be processed\r
+        /// \param gradeMV - the considered grade\r
+        /// \param moreThanOne - true if it the first element to display (should we put a '+' before)\r
+        template<typename U>\r
+        friend void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1>, unsigned int gradeMV, bool& moreThanOne);\r
+        /// \endcond // do not comment this functions\r
+\r
+/*\r
+        /// \cond DEV\r
+        /// \brief functions that enables to  to be able to extract a component of a multivector as a multivector. As an example, mv1.e(1,3) will create a multivector whose 1,3 component will be the component of mv1.\r
+        /// \tparam Arguments - denotes a variadic list of index\r
+        /// \param listIndices - the list of indices\r
+        /// \return a homogeneous multivector filled with zero except the listIndices index\r
+        /// \todo future work + handle the cases when the number of indices is higher than the dimension, and when the index is too high\r
+        template<typename... List>\r
+        Mvec e(List ...listIndices){\r
+            // operation on these arguments\r
+            // First determine the grade, simply use the variadic function sizeof...()\r
+            const int gd = sizeof...(listIndices); //\r
+\r
+            // Second, from these parameters, compute the index in the corresponding VectorXd\r
+            const int idx = (Binomial<algebraDimension,gd>::value-1) - computeIdxFromList(algebraDimension,gd,listIndices...);\r
+\r
+            Mvec res;\r
+            res.mvData[gd] = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(Binomial<algebraDimension,gd>::value);\r
+            res.mvData[gd](idx) = mvData[gd](idx);\r
+\r
+            return res;\r
+        }\r
+        /// \endcond // do not comment this functions\r
+*/\r
+\r
+\r
+        /// \brief functions that enables to extract one component of this multivector and initialize another mv with this component\r
+        /// \tparam Arguments - denotes a variadic list of index\r
+        /// \param grade : the considered grade\r
+        /// \param sizeOfKVector: C(dimension, grade)\r
+        /// \param indexInKvector: the considerd index\r
+        /// \return a new mv with the right component at the right place\r
+        Mvec extractOneComponent(const int grade, const int sizeOfKVector, const int indexInKvector) const {\r
+            Mvec mv;\r
+            auto itThisMV = this->findGrade(grade);\r
+            if(itThisMV==this->mvData.end()){\r
+                return mv;\r
+            }\r
+            mv = mv.componentToOne(grade,indexInKvector);\r
+            mv.mvData.begin()->vec.coeffRef(indexInKvector) = itThisMV->vec.coeff(indexInKvector);\r
+            return mv;\r
+        }\r
+\r
+\r
+        /// \brief set of functions that return multivectors containing only the value '1' for the specified component.\r
+project_multivector_one_component\r
+\r
+    };  // end of the class definition\r
+\r
+\r
+    /* ------------------------------------------------------------------------------------------------ */\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>::Mvec():gradeBitmap(0)\r
+    {}\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>::Mvec(const Mvec& mv) : mvData(mv.mvData), gradeBitmap(mv.gradeBitmap)\r
+    {\r
+        //std::cout << "copy constructor " << std::endl;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>::Mvec(Mvec<T>&& multivector) : mvData(std::move(multivector.mvData)), gradeBitmap(multivector.gradeBitmap)\r
+    {\r
+        // the move constructor\r
+        //std::cout << "move constructor" << std::endl;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename U>\r
+    Mvec<T>::Mvec(const Mvec<U> &mv) : gradeBitmap(mv.gradeBitmap)\r
+    {\r
+        for(auto it=mv.mvData.begin(); it != mv.mvData.end(); ++it){\r
+            Kvec<T> kvec;\r
+            kvec.grade = it->grade;\r
+            kvec.vec = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(binomialArray[it->grade]);\r
+            for(unsigned int i=0; i<it->vec.size(); ++i)\r
+                kvec.vec.coeffRef(i) = T(it->vec.coeff(i));\r
+            mvData.push_back(kvec);\r
+        }\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename U>\r
+    Mvec<T>::Mvec(const U val) {\r
+        if(val != U(0)) {\r
+            gradeBitmap = 1;\r
+            Kvec<T> kvec;\r
+            kvec.vec =Eigen::Matrix<T, Eigen::Dynamic, 1>(1);\r
+            kvec.grade =0;\r
+            mvData.push_back(kvec);\r
+            mvData.begin()->vec.coeffRef(0) = val;\r
+        }\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>::~Mvec()\r
+    {}\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>& Mvec<T>::operator=(const Mvec& mv){\r
+        if(&mv == this) return *this;\r
+        gradeBitmap = mv.gradeBitmap;\r
+        mvData = mv.mvData;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>& Mvec<T>::operator=(Mvec&& mv){\r
+        mvData = std::move(mv.mvData);\r
+        gradeBitmap = mv.gradeBitmap;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator+(const Mvec<T> &mv2) const {\r
+        Mvec<T> mv3(*this);\r
+        for(auto & itMv : mv2.mvData) {\r
+            auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);\r
+            it->vec += itMv.vec;\r
+        }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator+(const S &value, const Mvec<U> &mv){\r
+        return mv + value;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator+(const S &value) const {\r
+        Mvec<T> mv(*this);\r
+        if(value != T(0)) {\r
+            auto it = mv.createVectorXdIfDoesNotExist(0);\r
+            it->vec.coeffRef(0) += value;\r
+        }\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>& Mvec<T>::operator+=(const Mvec& mv){\r
+        *this = *this + mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> operator-(const Mvec<T> &mv) { // unary -\r
+        Mvec<T> mv2(mv);\r
+        for(auto & itMv : mv2.mvData)\r
+            itMv.vec = -itMv.vec;\r
+        return mv2;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator-(const Mvec<T> &mv2) const {\r
+        Mvec<T> mv3(*this);\r
+        for(auto & itMv : mv2.mvData) {\r
+            auto it = mv3.createVectorXdIfDoesNotExist(itMv.grade);\r
+            it->vec -= itMv.vec;\r
+        }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T>& Mvec<T>::operator-=(const Mvec& mv){\r
+        *this = *this - mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator-(const S &value, const Mvec<U> &mv){\r
+        return (-mv) + value;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator-(const S &value) const {\r
+        Mvec<T> mv(*this);\r
+        if(value != T(0)) {\r
+            auto it = mv.createVectorXdIfDoesNotExist(0);\r
+            it->vec.coeffRef(0) -= value;\r
+        }\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator^(const Mvec<T> &mv2) const {\r
+#if 0 // only recursive version\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // This version (with recursive call) is only faster than the standard recursive call if\r
+        // the multivector mv1 and mv2 are homogeneous or near from homogeneous.\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1\r
+            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2\r
+            {\r
+                unsigned int grade_mv3 = itMv1.grade + itMv2.grade;\r
+                if(grade_mv3 <= algebraDimension) {\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);\r
+                    outerProductHomogeneous(itMv1.vec, itMv2.vec, itMv3->vec,\r
+                                            itMv1.grade, itMv2.grade, grade_mv3);\r
+                }\r
+            }\r
+        return mv3;\r
+#else // use the adaptative pointer function array\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled outer function using the functions pointer called outerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+                if(itMv1.grade + itMv2.grade <= (int) algebraDimension ){\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(itMv1.grade + itMv2.grade);\r
+                    outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+                }\r
+            }\r
+        return mv3;\r
+#endif\r
+    }\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator^(const S &value, const Mvec<U> &mv) {\r
+        return  mv^value;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator^(const S &value) const {\r
+        Mvec<T> mv2(*this);\r
+        for(auto & itMv : mv2.mvData)\r
+            itMv.vec *= T(value);\r
+        return mv2;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> &Mvec<T>::operator^=(const Mvec &mv) {\r
+        *this = *this ^ mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator|(const Mvec<T> &mv2) const{\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+\r
+                // perform the inner product\r
+                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));\r
+                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);\r
+                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+\r
+                // check if the result is non-zero\r
+                if(!((itMv3->vec.array() != 0.0).any())){\r
+                    mv3.mvData.erase(itMv3);\r
+                    mv3.gradeBitmap &= ~(1<<absGradeMv3);\r
+                }\r
+            }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator|(const S &value, const Mvec<U> &mv) {\r
+        return mv | value;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator|(const S &value) const {\r
+        return (*this) > value ; // inner between a mv and a scalar gives 0 (empty multivector)\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> &Mvec<T>::operator|=(const Mvec &mv) {\r
+        *this = *this | mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator>(const Mvec<T> &mv2) const{\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+\r
+                // right contraction constraint: gradeMv1 >= gradeMv2\r
+                if(itMv1.grade < itMv2.grade)\r
+                    continue;\r
+\r
+                // perform the inner product\r
+                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));\r
+                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);\r
+                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+\r
+                // check if the result is non-zero\r
+                if(!((itMv3->vec.array() != 0.0).any())){\r
+                    mv3.mvData.erase(itMv3);\r
+                    mv3.gradeBitmap &= ~(1<<absGradeMv3);\r
+                }\r
+            }\r
+\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator>(const S &value, const Mvec<U> &mv) {\r
+        if( (mv.gradeBitmap & 1) == 0) return Mvec<U>();\r
+        else return Mvec<U>( U(mv.mvData.begin()->vec.coeff(0) * value ) );\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator>(const S &value) const {\r
+        // return mv x value\r
+        Mvec<T> mv(*this);\r
+        for(auto & itMv : mv.mvData)\r
+            itMv.vec *= T(value);\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator<(const Mvec<T> &mv2) const{\r
+\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+\r
+                // left contraction constraint: gradeMv1 <= gradeMv2\r
+                if(itMv1.grade > itMv2.grade)\r
+                    continue;\r
+\r
+                // perform the inner product\r
+                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));\r
+                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);\r
+                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+\r
+                // check if the result is non-zero\r
+                if(!((itMv3->vec.array() != 0.0).any())){\r
+                    mv3.mvData.erase(itMv3);\r
+                    mv3.gradeBitmap &= ~(1<<absGradeMv3);\r
+                }\r
+            }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator<(const S &value, const Mvec<U> &mv) {\r
+        // return mv x value\r
+        Mvec<U> mv3(mv);\r
+        for(auto & itMv : mv3.mvData)\r
+            itMv.vec *= U(value);\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator<(const S &value) const {\r
+        if( (gradeBitmap & 1) == 0) return Mvec<T>();\r
+        else return Mvec<T>( T(mvData.begin()->vec.coeff(0) * value ) );\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::scalarProduct(const Mvec<T> &mv2) const{\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+\r
+                if(itMv1.grade != itMv2.grade)\r
+                    continue;\r
+\r
+                // perform the inner product\r
+                int absGradeMv3 = 0;\r
+                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);\r
+                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+\r
+                // check if the result is non-zero\r
+                if(!((itMv3->vec.array() != 0.0).any())){\r
+                    mv3.mvData.erase(itMv3);\r
+                    mv3.gradeBitmap &= ~(1<<absGradeMv3);\r
+                }\r
+            }\r
+        return mv3;\r
+    }\r
+\r
+\r
+project_singular_metric_comment_begin\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::outerPrimalDual(const Mvec<T> &mv2) const{\r
+        Mvec<T> mv3;\r
+\r
+        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1\r
+            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2\r
+            {\r
+                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);\r
+                if(grade_mv3 <= algebraDimension) {\r
+                    // handle the scalar product as well as the left contraction\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);\r
+                    outerProductPrimalDual(itMv1.vec, itMv2.vec, itMv3->vec,\r
+                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));\r
+\r
+                    // check if the result is non-zero\r
+                    if(!((itMv3->vec.array() != 0.0).any())){\r
+                        mv3.mvData.erase(itMv3);\r
+                        mv3.gradeBitmap &= ~(1<<grade_mv3);\r
+                    }\r
+                }\r
+            }\r
+\r
+        return mv3;\r
+    }\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::outerDualPrimal(const Mvec<T> &mv2) const{\r
+        Mvec<T> mv3;\r
+\r
+        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1\r
+            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2\r
+            {\r
+                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);\r
+                if(grade_mv3 <= algebraDimension) {\r
+                    // handle the scalar product as well as the left contraction\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);\r
+                    outerProductDualPrimal(itMv1.vec, itMv2.vec, itMv3->vec,\r
+                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));\r
+\r
+                    // check if the result is non-zero\r
+                    if(!((itMv3->vec.array() != 0.0).any())){\r
+                        mv3.mvData.erase(itMv3);\r
+                        mv3.gradeBitmap &= ~(1<<grade_mv3);\r
+                    }\r
+                }\r
+            }\r
+\r
+\r
+        return mv3;\r
+\r
+    }\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::outerDualDual(const Mvec<T> &mv2) const{\r
+        Mvec<T> mv3;\r
+\r
+        for(const auto & itMv1 : this->mvData)    // all per-grade component of mv1\r
+            for(const auto & itMv2 : mv2.mvData)  // all per-grade component of mv2\r
+            {\r
+                unsigned int grade_mv3 = itMv1.grade + (algebraDimension-itMv2.grade);\r
+                if(grade_mv3 <= algebraDimension) {\r
+                    // handle the scalar product as well as the left contraction\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(grade_mv3);\r
+                    outerProductDualDual(itMv1.vec, itMv2.vec, itMv3->vec,\r
+                                           itMv1.grade, itMv2.grade, (unsigned)(algebraDimension-grade_mv3));\r
+\r
+                    // check if the result is non-zero\r
+                    if(!((itMv3->vec.array() != 0.0).any())){\r
+                        mv3.mvData.erase(itMv3);\r
+                        mv3.gradeBitmap &= ~(1<<grade_mv3);\r
+                    }\r
+                }\r
+            }\r
+\r
+\r
+        return mv3;\r
+\r
+    }\r
+\r
+project_singular_metric_comment_end\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::hestenesProduct(const Mvec<T> &mv2) const{\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled inner function using the functions pointer called innerFunctionsContainer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+\r
+                if(itMv1.grade*itMv2.grade == 0)\r
+                    continue;\r
+\r
+                // perform the inner product\r
+                int absGradeMv3 = std::abs((int)(itMv1.grade - itMv2.grade));\r
+                auto itMv3 = mv3.createVectorXdIfDoesNotExist(absGradeMv3);\r
+                innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+\r
+                // check if the result is non-zero\r
+                if(!((itMv3->vec.array() != 0.0).any())){\r
+                    mv3.mvData.erase(itMv3);\r
+                    mv3.gradeBitmap &= ~(1<<absGradeMv3);\r
+                }\r
+            }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator*(const Mvec<T> &mv2) const {\r
+        // Loop over non-empty grade of mv1 and mv2\r
+        // call the right explicit unrolled product function using the functions pointer\r
+        Mvec<T> mv3;\r
+        for(const auto & itMv1 : this->mvData)\r
+            for(const auto & itMv2 : mv2.mvData){\r
+                project_select_recursive_geometric_product_template\r
+                // outer product block\r
+                unsigned int gradeOuter = itMv1.grade + itMv2.grade;\r
+                if(gradeOuter <=  algebraDimension ){\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeOuter);\r
+                    outerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+                    if(!((itMv3->vec.array() != 0.0).any())){\r
+                        mv3.mvData.erase(itMv3);\r
+                        mv3.gradeBitmap &= ~(1<<gradeOuter);\r
+                    }\r
+                }\r
+\r
+                // inner product block\r
+                unsigned int gradeInner = (unsigned int)std::abs((int)(itMv1.grade-itMv2.grade));\r
+                // when the grade of one of the kvectors is zero, the inner product is the same as the outer product\r
+                if(gradeInner != gradeOuter) {\r
+                    auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeInner);\r
+                    innerFunctionsContainer<T>[itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+                    // check if the result is non-zero\r
+                    if(!((itMv3->vec.array() != 0.0).any())){\r
+                        mv3.mvData.erase(itMv3);\r
+                        mv3.gradeBitmap &= ~(1<<gradeInner);\r
+                    }\r
+\r
+                    // geometric product part\r
+                    int gradeMax = std::min(((2*algebraDimension)-gradeOuter)+1,gradeOuter);\r
+                    for (int gradeResult = gradeInner+2; gradeResult < gradeMax; gradeResult+=2) {\r
+                        auto itMv3 = mv3.createVectorXdIfDoesNotExist(gradeResult);\r
+                        geometricFunctionsContainer<T>[gradeResult][itMv1.grade][itMv2.grade](itMv1.vec, itMv2.vec, itMv3->vec);\r
+                        // check if the result is non-zero\r
+                        if(!((itMv3->vec.array() != 0.0).any())){\r
+                            mv3.mvData.erase(itMv3);\r
+                            mv3.gradeBitmap &= ~(1<<gradeResult);\r
+                        }\r
+                    }\r
+                }\r
+            }\r
+        return mv3;\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator*(const S &value, const Mvec<U> &mv){\r
+        return mv * value;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator*(const S &value) const{\r
+        Mvec<T> mv(*this);\r
+        for(auto & itMv : mv.mvData)\r
+            itMv.vec *= T(value);\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> &Mvec<T>::operator*=(const Mvec &mv) {\r
+        *this = *this * mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::operator/(const Mvec<T> &mv2) const {\r
+        return *this * mv2.inv();\r
+    }\r
+\r
+\r
+    template<typename U, typename S>\r
+    Mvec<U> operator/(const S &value, const Mvec<U> &mv){\r
+        return value * mv.inv();\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    template<typename S>\r
+    Mvec<T> Mvec<T>::operator/(const S &value) const {\r
+        Mvec<T> mv(*this);\r
+        for(auto & itMv : mv.mvData)\r
+            itMv.vec /= value;\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> &Mvec<T>::operator/=(const Mvec &mv) {\r
+        *this = *this / mv;\r
+        return *this;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> operator~(const Mvec<T> &mv){\r
+        return mv.reverse();\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::inv() const {\r
+        T n = this->quadraticNorm();\r
+        if(n<std::numeric_limits<T>::epsilon() && n>-std::numeric_limits<T>::epsilon())\r
+            return Mvec<T>(); // return 0, this is was gaviewer does.\r
+        return this->reverse() / n;\r
+    };\r
+\r
+\r
+    project_singular_metric_comment_begin\r
+    template<typename U>\r
+    Mvec<U> operator!(const Mvec<U> &mv) {\r
+        return mv.dual();\r
+    }\r
+project_singular_metric_comment_end\r
+\r
+    /// \brief defines the left contraction between two multivectors\r
+    /// \param mv1 - a multivector\r
+    /// \param mv2 - second operand, corresponds to a multivector\r
+    /// \return mv1 left contraction mv2\r
+    template<typename T>\r
+    Mvec<T> leftContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){\r
+        return mv1 < mv2;\r
+    }\r
+\r
+\r
+    /// \brief defines the right contraction between two multivectors\r
+    /// \param mv1 - a multivector\r
+    /// \param mv2 - second operand, corresponds to a multivector\r
+    /// \return mv1 right contraction mv2\r
+    template<typename T>\r
+    Mvec<T> rightContraction(const Mvec<T> &mv1, const Mvec<T> &mv2){\r
+        return mv1 > mv2;\r
+    }\r
+\r
+\r
+    /// \brief returns a multivector that only contains the coefficient associated to the pseudoscalar.\r
+    /// \return an empty Mvec if the requested element is not part of the multivector, or the multivector that contains only this element if present in the current multivector.\r
+    template<typename T>\r
+    Mvec<T> I() {\r
+        Mvec<T> mvec;\r
+        mvec[project_pseudo_scalar] = 1.0;\r
+        return mvec;\r
+    }\r
+\r
+project_singular_metric_comment_begin\r
+    /// \brief return the inverse of the pseudo scalar.\r
+    /// \return returns a multivector corresponding to the inverse of the pseudo scalar.\r
+    template<typename T>\r
+    Mvec<T> Iinv() {\r
+        Mvec<T> mvec;\r
+        mvec[project_pseudo_scalar] = pseudoScalarInverse; // we do not return just a scalar of type T because the grade should be dimension and not 0.\r
+        return mvec;\r
+    }\r
+project_singular_metric_comment_end\r
+\r
+project_singular_metric_comment_begin\r
+    // compute the dual of a multivector (i.e mv* = mv.reverse() * Iinv)\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::dual() const {\r
+        Mvec<T> mvResult;\r
+        // for each k-vectors of the multivector\r
+        for(auto itMv=mvData.rbegin(); itMv!=mvData.rend(); ++itMv){\r
+\r
+            // create the dual k-vector\r
+            Kvec<T> kvec ={itMv->vec,algebraDimension-(itMv->grade)};\r
+\r
+            // some elements need to be permuted\r
+            for(unsigned int i=0;i<binomialArray[itMv->grade];++i)\r
+                kvec.vec.coeffRef(dualPermutations[itMv->grade][i]) = itMv->vec.coeff(i);\r
+\r
+            // the inner product may involve some constant multiplucation for the dual elements\r
+            kvec.vec = kvec.vec.cwiseProduct(dualCoefficients[itMv->grade].template cast<T>());\r
+\r
+            // add the k-vector to the resulting multivector\r
+            mvResult.mvData.push_back(kvec);\r
+            mvResult.gradeBitmap |= (1 << kvec.grade);\r
+        }\r
+        return mvResult;\r
+    };\r
+project_singular_metric_comment_end\r
+\r
+    // \brief compute the reverse of a multivector\r
+    // \return - the reverse of the multivector\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::reverse() const {\r
+        Mvec<T> mv(*this);\r
+        for(auto & itMv : mv.mvData)\r
+            if(signReversePerGrade[itMv.grade] == -1)\r
+                itMv.vec *= -1;\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    Mvec<T> Mvec<T>::grade(const int i) const{\r
+\r
+        auto it = findGrade(i);\r
+        Mvec<T> mv;\r
+\r
+        // if not found, return the empty multivector\r
+        if(it == mvData.end())\r
+            return mv;\r
+\r
+        // else return the grade 'i' data\r
+        mv.mvData.push_back(*it);\r
+        mv.gradeBitmap = 1 << (i);\r
+\r
+        return mv;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    std::vector<unsigned int> Mvec<T>::grades() const {\r
+        std::vector<unsigned int> mvGrades;\r
+        for(unsigned int i=0; i< algebraDimension+1; ++i){\r
+\r
+            if( (gradeBitmap & (1<<i)))\r
+                mvGrades.push_back(i);\r
+        }\r
+\r
+        return mvGrades;\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    void Mvec<T>::roundZero(const T epsilon) {\r
+        // loop over each k-vector of the multivector\r
+        auto itMv=mvData.begin();\r
+        while(itMv != mvData.end()){\r
+            // loop over each element of the k-vector\r
+            for(unsigned int i=0; i<(unsigned int)itMv->vec.size(); ++i)\r
+                if(fabs(itMv->vec.coeff(i)) <= epsilon)\r
+                    itMv->vec.coeffRef(i) = 0.0;\r
+            // if the k-vector is full of 0, remove it\r
+            if(!((itMv->vec.array() != 0.0).any())){\r
+                gradeBitmap = gradeBitmap - (1 << itMv->grade);\r
+                mvData.erase(itMv++);\r
+            }\r
+            else ++itMv;\r
+        }\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    void Mvec<T>::clear(const int grade) {\r
+\r
+        // full erase\r
+        if(grade < 0){\r
+            mvData.clear();\r
+            gradeBitmap = 0;\r
+            return;\r
+        }\r
+\r
+        // partial erase\r
+        auto iter = findGrade(grade);\r
+        if(iter != mvData.end()) {\r
+            gradeBitmap = gradeBitmap - (1 << iter->grade);\r
+            iter = mvData.erase(iter);\r
+        }\r
+    }\r
+\r
+\r
+    template<typename T>\r
+    void Mvec<T>::display() const {\r
+        if(mvData.size() == 0)\r
+            std::cout << " null grade , null value " <<std::endl;\r
+\r
+        for(auto itmv=mvData.begin(); itmv!=mvData.end(); ++itmv){\r
+            std::cout << "  grade   : " << itmv->grade  << std::endl;\r
+            std::cout << "  kvector : " << itmv->vec.transpose() << std::endl;\r
+        }\r
+        std::cout << std::endl;\r
+    }\r
+\r
+\r
+    /// \cond DEV\r
+    template<typename U>\r
+    void traverseKVector(std::ostream &stream, const Eigen::Matrix<U, Eigen::Dynamic, 1> kvector, unsigned int gradeMV, bool& moreThanOne ){\r
+\r
+        // version with XOR indices\r
+        // for the current grade, generate all the XOR indices\r
+        std::vector<bool> booleanXorIndex(algebraDimension);\r
+        std::fill(booleanXorIndex.begin(), booleanXorIndex.end(), false);\r
+        // build the first xorIndex of grade 'grade' (with 'grade' values to true).\r
+        std::fill(booleanXorIndex.begin(), booleanXorIndex.begin() + gradeMV, true);\r
+        unsigned int positionInKVector = 0;\r
+\r
+        do {\r
+            // convert the vector of bool to the string containing all the basis vectors\r
+            std::string basisBlades={};\r
+            for(unsigned int i=0; i<algebraDimension; ++i) {\r
+                if(booleanXorIndex[i]) {\r
+                    basisBlades += basisVectors[i];\r
+                }\r
+            }\r
+\r
+            if(kvector.coeff(positionInKVector)!= 0) {\r
+                if(!(moreThanOne)){\r
+                    stream<< kvector.coeff(positionInKVector) << "*e"+ basisBlades;\r
+                    moreThanOne = true;\r
+                }else{\r
+                    if(kvector.coeff(positionInKVector)>0)\r
+                        stream<< " + " << kvector.coeff(positionInKVector) << "*e" + basisBlades;\r
+                    else\r
+                        stream<< " - " << -kvector.coeff(positionInKVector) << "*e" + basisBlades;\r
+                }\r
+            }\r
+\r
+            positionInKVector++;\r
+\r
+        } while(std::prev_permutation(booleanXorIndex.begin(), booleanXorIndex.end())); // compute next permutation of the true values of booleanXorIndex\r
+\r
+    }\r
+    /// \endcond\r
+\r
+\r
+    template<typename U>\r
+    std::ostream &operator<<(std::ostream &stream, const Mvec<U> &mvec) {\r
+        if(mvec.mvData.size()==0){\r
+            stream << " 0 ";\r
+            return stream;\r
+        }\r
+\r
+        bool moreThanOne = false;\r
+        auto mvIterator = mvec.mvData.begin();\r
+\r
+        // if the iterator corresponds to the scalar\r
+        if(mvIterator->grade == 0){\r
+            stream << mvIterator->vec.coeff(0);\r
+            moreThanOne = true;\r
+            ++mvIterator;\r
+        }\r
+\r
+        // for all other k-vectors of mvec\r
+        for(; mvIterator!=mvec.mvData.end(); ++mvIterator){\r
+\r
+            // call the function that covers a single k-vector\r
+            traverseKVector(stream,mvIterator->vec,mvIterator->grade, moreThanOne);\r
+        }\r
+\r
+        if(!moreThanOne)\r
+            stream << " 0 ";\r
+\r
+        return stream;\r
+    }\r
+\r
+    // static unit multivector functions (returns a multivector with only one component)\r
+project_static_multivector_one_component\r
+\r
+    //    template<typename U>\r
+    //    void recursiveTraversalMultivector(std::ostream &stream, const Mvec<U> &mvec, unsigned int currentGrade, int currentIndex, std::vector<int> listBasisBlades, unsigned int lastIndex, unsigned int gradeMV, bool& moreThanOne);\r
+\r
+\r
+    void temporaryFunction1();\r
+\r
+}     /// End of Namespace\r
+\r
+#endif // project_inclusion_guard\r
index 5100fc3..f13e501 100644 (file)
@@ -1,39 +1,39 @@
-#include <iostream>
-#include <cstdlib>
-#include <project_namespace/Mvec.hpp>
-
-
-
-int main(){
-
-    // sample instructions
-    std::cout << "metric : \n" << project_namespace::metric << std::endl;
-
-    // accessor
-    project_namespace::Mvec<double> mv1;
-    mv1[project_namespace::scalar] = 1.0;
-    mv1[project_namespace::Eproject_first_vector_basis] = 42.0;
-    std::cout << "mv1 : " << mv1 << std::endl;
-
-    project_namespace::Mvec<double> mv2;
-    mv2[project_namespace::Eproject_first_vector_basis] = 1.0;
-    mv2[project_namespace::Eproject_second_vector_basis] = 2.0;
-    mv2 += project_namespace::I<double>() + 2*project_namespace::eproject_first_vector_basisproject_second_vector_basis<double>();
-    std::cout << "mv2 : " << mv2 << std::endl << std::endl;
-
-    // some products
-    std::cout << "outer product     : " << (mv1 ^ mv2) << std::endl;
-    std::cout << "inner product     : " << (mv1 | mv2) << std::endl;
-    std::cout << "geometric product : " << (mv1 * mv2) << std::endl;
-    std::cout << "left contractiont : " << (mv1 < mv2) << std::endl;
-    std::cout << "right contractiont: " << (mv1 > mv2) << std::endl;
-    std::cout << std::endl;
-
-    // some tools
-    std::cout << "grade : " << mv1.grade()  << std::endl;
-    std::cout << "norm  : " << mv1.norm()  << std::endl;
-    mv1.clear();
-    if(mv1.isEmpty()) std::cout << "mv1 is empty: ok" << std::endl;
-
-    return EXIT_SUCCESS;
-}
+#include <iostream>\r
+#include <cstdlib>\r
+#include <project_namespace/Mvec.hpp>\r
+\r
+\r
+\r
+int main(){\r
+\r
+    // sample instructions\r
+    std::cout << "metric : \n" << project_namespace::metric << std::endl;\r
+\r
+    // accessor\r
+    project_namespace::Mvec<double> mv1;\r
+    mv1[project_namespace::scalar] = 1.0;\r
+    mv1[project_namespace::Eproject_first_vector_basis] = 42.0;\r
+    std::cout << "mv1 : " << mv1 << std::endl;\r
+\r
+    project_namespace::Mvec<double> mv2;\r
+    mv2[project_namespace::Eproject_first_vector_basis] = 1.0;\r
+    mv2[project_namespace::Eproject_second_vector_basis] = 2.0;\r
+    mv2 += project_namespace::I<double>() + 2*project_namespace::eproject_first_vector_basisproject_second_vector_basis<double>();\r
+    std::cout << "mv2 : " << mv2 << std::endl << std::endl;\r
+\r
+    // some products\r
+    std::cout << "outer product     : " << (mv1 ^ mv2) << std::endl;\r
+    std::cout << "inner product     : " << (mv1 | mv2) << std::endl;\r
+    std::cout << "geometric product : " << (mv1 * mv2) << std::endl;\r
+    std::cout << "left contraction : " << (mv1 < mv2) << std::endl;\r
+    std::cout << "right contraction: " << (mv1 > mv2) << std::endl;\r
+    std::cout << std::endl;\r
+\r
+    // some tools\r
+    std::cout << "grade : " << mv1.grade()  << std::endl;\r
+    std::cout << "norm  : " << mv1.norm()  << std::endl;\r
+    mv1.clear();\r
+    if(mv1.isEmpty()) std::cout << "mv1 is empty: ok" << std::endl;\r
+\r
+    return EXIT_SUCCESS;\r
+}\r
index 7a374d9..7dff4eb 100644 (file)
-// Copyright (c) 2018 by University Paris-Est Marne-la-Vallee
-// Directory.cpp
-// This file is part of the Garamon Generator.
-// Authors: Stephane Breuils and Vincent Nozick
-// Contact: vincent.nozick@u-pem.fr
-//
-// Licence MIT
-// A a copy of the MIT License is given along with this program
-
-
-#include "Directory.hpp"
-
-#include <iostream>
-#include <cstdlib>
-#include <fstream>
-#include <sstream>
-#include <sys/types.h> // required for stat.h
-#include <sys/stat.h>
-#include <regex>       // for substitute
-
-
-#if defined(_WIN32) && defined(__MINGW32__)
-       #include <direct.h>
-       #include <io.h>
-#endif // __WINDOWS Mingw compiler__
-
-
-void makeDirectory(const std::string &dirName) {
-
-    int nError = 0;
-
-#if defined(_WIN32)
-       nError = _mkdir(dirName.c_str()); // can be used on Windows
-#else
-    mode_t nMode = 0733; // UNIX style permissions
-    nError = mkdir(dirName.c_str(), nMode); // can be used on non-Windows
-#endif
-
-    // handle your error here
-    if (nError != 0) {
-        std::cerr << "error: can not create directory: " << dirName << std::endl;
-        exit(EXIT_FAILURE);
-    }
-}
-
-
-bool directoryExists(const std::string &dirName) {
-
-    struct stat info;
-
-    if(stat( dirName.c_str(), &info ) != 0)
-        return false; // cannot access to the directory
-
-    if(info.st_mode & S_IFDIR) // in case of problem, check with: S_ISDIR()
-        return true;
-
-    return false;  // file exists but is not a directory
-}
-
-
-bool directoryOrFileExists(const std::string &dirName) {
-
-    struct stat info;
-
-    if(stat( dirName.c_str(), &info ) == 0)
-        return true; // can access to the directory or file
-
-    return false;
-}
-
-
-bool directoryOrFileExists_ifstream (const std::string& name) {
-    std::ifstream f(name.c_str());
-    return f.good();
-}
-
-
-std::string readFile(const std::string &fileName) {
-
-    std::ifstream myfile;
-    myfile.open(fileName, std::ios::in);
-
-    // check if the file is opened
-    if(!myfile.is_open()){
-        std::cerr << "error: can not open file: " << fileName << std::endl;
-        exit(EXIT_FAILURE);
-    }
-
-    // copy the data to a string
-    std::string data;
-    data.assign( (std::istreambuf_iterator<char>(myfile)) , (std::istreambuf_iterator<char>()) );
-
-    myfile.close();
-
-    return data;
-}
-
-
-bool writeFile(const std::string &data, const std::string &fileName) {
-
-    std::ofstream myfile(fileName);
-    if(!myfile.is_open()){
-        std::cerr << "error: can not create file: " << fileName << std::endl;
-        return false;
-    }
-
-    myfile << data;
-    myfile.close();
-
-    return true;
-}
-
-
-
-
-
-void substitute(std::string &data, const std::string &pattern, const std::string &replaceBy) {
-
-    data = std::regex_replace(data, std::regex(pattern), replaceBy);
-}
-
-
-bool copyBin(const std::string &src, const std::string &dest) {
-
-    std::ifstream srcFile(src, std::ios::binary);
-    std::ofstream destFile(dest, std::ios::binary);
-    destFile << srcFile.rdbuf();
-
-    return srcFile && destFile;
-}
-
-
-bool copyText(const std::string &src, const std::string &dest) {
-
-    std::ifstream srcFile(src);
-    std::ofstream destFile(dest);
-    destFile << srcFile.rdbuf();
-
-    return srcFile && destFile;
-}
+// Copyright (c) 2018 by University Paris-Est Marne-la-Vallee\r
+// Directory.cpp\r
+// This file is part of the Garamon Generator.\r
+// Authors: Stephane Breuils and Vincent Nozick\r
+// Contact: vincent.nozick@u-pem.fr\r
+//\r
+// Licence MIT\r
+// A a copy of the MIT License is given along with this program\r
+\r
+\r
+#include "Directory.hpp"\r
+\r
+#include <iostream>\r
+#include <cstdlib>\r
+#include <fstream>\r
+#include <sstream>\r
+#include <sys/types.h> // required for stat.h\r
+#include <sys/stat.h>\r
+#include <regex>       // for substitute\r
+\r
+\r
+#if defined(_WIN32) && (defined(__MINGW32__) || defined(_MSC_BUILD))\r
+       #include <direct.h>\r
+       #include <io.h>\r
+#endif // __WINDOWS Mingw or MSVC compiler__\r
+\r
+\r
+void makeDirectory(const std::string &dirName) {\r
+\r
+    int nError = 0;\r
+\r
+#if defined(_WIN32)\r
+       nError = _mkdir(dirName.c_str()); // can be used on Windows\r
+#else\r
+    mode_t nMode = 0733; // UNIX style permissions\r
+    nError = mkdir(dirName.c_str(), nMode); // can be used on non-Windows\r
+#endif\r
+\r
+    // handle your error here\r
+    if (nError != 0) {\r
+        std::cerr << "error: can not create directory: " << dirName << std::endl;\r
+        exit(EXIT_FAILURE);\r
+    }\r
+}\r
+\r
+\r
+bool directoryExists(const std::string &dirName) {\r
+\r
+    struct stat info;\r
+\r
+    if(stat( dirName.c_str(), &info ) != 0)\r
+        return false; // cannot access to the directory\r
+\r
+    if(info.st_mode & S_IFDIR) // in case of problem, check with: S_ISDIR()\r
+        return true;\r
+\r
+    return false;  // file exists but is not a directory\r
+}\r
+\r
+\r
+bool directoryOrFileExists(const std::string &dirName) {\r
+\r
+    struct stat info;\r
+\r
+    if(stat( dirName.c_str(), &info ) == 0)\r
+        return true; // can access to the directory or file\r
+\r
+    return false;\r
+}\r
+\r
+\r
+bool directoryOrFileExists_ifstream (const std::string& name) {\r
+    std::ifstream f(name.c_str());\r
+    return f.good();\r
+}\r
+\r
+\r
+std::string readFile(const std::string &fileName) {\r
+\r
+    std::ifstream myfile;\r
+    myfile.open(fileName, std::ios::in);\r
+\r
+    // check if the file is opened\r
+    if(!myfile.is_open()){\r
+        std::cerr << "error: can not open file: " << fileName << std::endl;\r
+        exit(EXIT_FAILURE);\r
+    }\r
+\r
+    // copy the data to a string\r
+    std::string data;\r
+    data.assign( (std::istreambuf_iterator<char>(myfile)) , (std::istreambuf_iterator<char>()) );\r
+\r
+    myfile.close();\r
+\r
+    return data;\r
+}\r
+\r
+\r
+bool writeFile(const std::string &data, const std::string &fileName) {\r
+\r
+    std::ofstream myfile(fileName);\r
+    if(!myfile.is_open()){\r
+        std::cerr << "error: can not create file: " << fileName << std::endl;\r
+        return false;\r
+    }\r
+\r
+    myfile << data;\r
+    myfile.close();\r
+\r
+    return true;\r
+}\r
+\r
+\r
+\r
+\r
+\r
+void substitute(std::string &data, const std::string &pattern, const std::string &replaceBy) {\r
+\r
+    data = std::regex_replace(data, std::regex(pattern), replaceBy);\r
+}\r
+\r
+\r
+bool copyBin(const std::string &src, const std::string &dest) {\r
+\r
+    std::ifstream srcFile(src, std::ios::binary);\r
+    std::ofstream destFile(dest, std::ios::binary);\r
+    destFile << srcFile.rdbuf();\r
+\r
+    return srcFile && destFile;\r
+}\r
+\r
+\r
+bool copyText(const std::string &src, const std::string &dest) {\r
+\r
+    std::ifstream srcFile(src);\r
+    std::ofstream destFile(dest);\r
+    destFile << srcFile.rdbuf();\r
+\r
+    return srcFile && destFile;\r
+}\r
index 7b1ac89..21f1549 100644 (file)
 #include <Eigen/Dense> // use the metric which is dense\r
 #include <Eigen/Sparse> // use the transformation matrices which are sparse\r
 #include <iostream>\r
+#if defined(_MSC_BUILD)\r
+       #ifndef __builtin_popcount\r
+               #define __builtin_popcount __popcnt\r
+       #endif\r
+#endif // __WINDOWS MSVC compiler__\r
 \r
 \r
 /// In the computation of mv3 = mv1^mv2, this represents a quadruplet containing the indices of mv1,mv2,mv3 and a coefficient in a product between two blades.\r