/Users/barrand/private/dev/softinex/old/inexlib-1.2/inlib/inlib/a4

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file inlib.license for terms.

#ifndef inlib_a4
#define inlib_a4

#include "sqm"

#include <cmath> //sqrt

namespace inlib {
namespace a4 {

// vec is an a4

template <class T>
class vec : public array<T> {
public:
  vec():array<T>(1,4){}
  vec(const T a_vec[4]):array<T>(1,4) {
    array<T>::m_vector[0] = a_vec[0];
    array<T>::m_vector[1] = a_vec[1];
    array<T>::m_vector[2] = a_vec[2];
    array<T>::m_vector[3] = a_vec[3];
  }
  vec(const T& a0,const T& a1,const T& a2,const T& a3):array<T>(1,4) {
    array<T>::m_vector[0] = a0;
    array<T>::m_vector[1] = a1;
    array<T>::m_vector[2] = a2;
    array<T>::m_vector[3] = a3;
  }
  virtual ~vec() {}
public:
  vec(const vec& a_from):array<T>(a_from){}
  vec& operator=(const vec& a_from) {
    array<T>::operator=(a_from);
    return *this;
  }
public:
  void set_value(const T& a1,const T& a2,const T& a3,const T& a4) {
    array<T>::m_vector[0] = a1;
    array<T>::m_vector[1] = a2;
    array<T>::m_vector[2] = a3;
    array<T>::m_vector[3] = a4;
  }
  T v0() const { return array<T>::m_vector[0];}
  T v1() const { return array<T>::m_vector[1];}
  T v2() const { return array<T>::m_vector[2];}
  T v3() const { return array<T>::m_vector[3];}

  void v0(const T& a_value) { array<T>::m_vector[0] = a_value;}
  void v1(const T& a_value) { array<T>::m_vector[1] = a_value;}
  void v2(const T& a_value) { array<T>::m_vector[2] = a_value;}
  void v3(const T& a_value) { array<T>::m_vector[3] = a_value;}

  bool set_value(unsigned int a_index,const T& a_value) { 
    if(a_index>=4) return false;
    array<T>::m_vector[a_index] = a_value;
    return true;
  }

  T length() const {
    T m0 = array<T>::m_vector[0];
    T m1 = array<T>::m_vector[1];
    T m2 = array<T>::m_vector[2];
    T m3 = array<T>::m_vector[3];
    return (T)::sqrt(m0*m0+m1*m1+m2*m2+m3*m3);
  }
  T normalize() {
    T norme = length();
    if(norme==T()) return 0;
    divide(norme);
    return norme;
  }
};

// sqm is an a44

template <class T>
class sqm : public inlib::sqm<T> {
public:
  sqm(): inlib::sqm<T>(4){}
  virtual ~sqm() {}
public:
  sqm(const sqm& a_from): inlib::sqm<T>(a_from){}
  sqm& operator=(const sqm& a_from) {
    inlib::sqm<T>::operator=(a_from);
    return *this;
  }
public:
  void v00(const T& a_value){array<T>::m_vector[0+0*4] = a_value;}
  void v10(const T& a_value){array<T>::m_vector[1+0*4] = a_value;}
  void v20(const T& a_value){array<T>::m_vector[2+0*4] = a_value;}
  void v30(const T& a_value){array<T>::m_vector[3+0*4] = a_value;}

  void v01(const T& a_value){array<T>::m_vector[0+1*4] = a_value;}
  void v11(const T& a_value){array<T>::m_vector[1+1*4] = a_value;}
  void v21(const T& a_value){array<T>::m_vector[2+1*4] = a_value;}
  void v31(const T& a_value){array<T>::m_vector[3+1*4] = a_value;}

  void v02(const T& a_value){array<T>::m_vector[0+2*4] = a_value;}
  void v12(const T& a_value){array<T>::m_vector[1+2*4] = a_value;}
  void v22(const T& a_value){array<T>::m_vector[2+2*4] = a_value;}
  void v32(const T& a_value){array<T>::m_vector[3+2*4] = a_value;}

  void v03(const T& a_value){array<T>::m_vector[0+3*4] = a_value;}
  void v13(const T& a_value){array<T>::m_vector[1+3*4] = a_value;}
  void v23(const T& a_value){array<T>::m_vector[2+3*4] = a_value;}
  void v33(const T& a_value){array<T>::m_vector[3+3*4] = a_value;}

  T v00() const {return array<T>::m_vector[0+0*4];}
  T v10() const {return array<T>::m_vector[1+0*4];}
  T v20() const {return array<T>::m_vector[2+0*4];}
  T v30() const {return array<T>::m_vector[3+0*4];}

  T v01() const {return array<T>::m_vector[0+1*4];}
  T v11() const {return array<T>::m_vector[1+1*4];}
  T v21() const {return array<T>::m_vector[2+1*4];}
  T v31() const {return array<T>::m_vector[3+1*4];}

  T v02() const {return array<T>::m_vector[0+2*4];}
  T v12() const {return array<T>::m_vector[1+2*4];}
  T v22() const {return array<T>::m_vector[2+2*4];}
  T v32() const {return array<T>::m_vector[3+2*4];}

  T v03() const {return array<T>::m_vector[0+3*4];}
  T v13() const {return array<T>::m_vector[1+3*4];}
  T v23() const {return array<T>::m_vector[2+3*4];}
  T v33() const {return array<T>::m_vector[3+3*4];}

  virtual T determinant() const {
    T e00 =  D33(v11(),v12(),v13(),
                      v21(),v22(),v23(),
                      v31(),v32(),v33());
    T e01 = -D33(v01(),v02(),v03(),
                      v21(),v22(),v23(),
                      v31(),v32(),v33());
    T e02 =  D33(v01(),v02(),v03(),
                      v11(),v12(),v13(),
                      v31(),v32(),v33());
    T e03 = -D33(v01(),v02(),v03(),
                      v11(),v12(),v13(),
                      v21(),v22(),v23());
                  
    return (  v00() * e00
            + v10() * e01
            + v20() * e02
            + v30() * e03);
  }

  bool invert(sqm<T>& a_result) const {
    // Twice faster than the generic invertion method.
    a_result.v00(
      D33(v11(),v12(),v13(),
          v21(),v22(),v23(),
        v31(),v32(),v33()));
    a_result.v01(
      -D33(v01(),v02(),v03(),
         v21(),v22(),v23(),
           v31(),v32(),v33()));
    a_result.v02(
      D33(v01(),v02(),v03(),
          v11(),v12(),v13(),
        v31(),v32(),v33()));
    a_result.v03(
     -D33(v01(),v02(),v03(),
          v11(),v12(),v13(),
        v21(),v22(),v23()));
                    
    T deter =
        v00() * a_result.v00()
      + v10() * a_result.v01()
      + v20() * a_result.v02()
      + v30() * a_result.v03();
  
    if(deter==array<T>::zero()) return false;
                    
    a_result.v10(
     -D33(v10(),v12(),v13(),
          v20(),v22(),v23(),
          v30(),v32(),v33()));
    a_result.v11(
      D33(v00(),v02(),v03(),
        v20(),v22(),v23(),
        v30(),v32(),v33()));
    a_result.v12(
     -D33(v00(),v02(),v03(),
          v10(),v12(),v13(),
        v30(),v32(),v33()));
    a_result.v13(
      D33(v00(),v02(),v03(),
        v10(),v12(),v13(),
        v20(),v22(),v23()));
                    
    a_result.v20(
      D33(v10(),v11(),v13(),
        v20(),v21(),v23(),
        v30(),v31(),v33()));
    a_result.v21(
     -D33(v00(),v01(),v03(),
        v20(),v21(),v23(),
        v30(),v31(),v33()));
    a_result.v22(
      D33(v00(),v01(),v03(),
        v10(),v11(),v13(),
        v30(),v31(),v33()));
    a_result.v23(
     -D33(v00(),v01(),v03(),
        v10(),v11(),v13(),
        v20(),v21(),v23()));
                    
    a_result.v30(
     -D33(v10(),v11(),v12(),
        v20(),v21(),v22(),
        v30(),v31(),v32()));
    a_result.v31(
      D33(v00(),v01(),v02(),
        v20(),v21(),v22(),
        v30(),v31(),v32()));
    a_result.v32(
     -D33(v00(),v01(),v02(),
        v10(),v11(),v12(),
        v30(),v31(),v32()));
    a_result.v33(
      D33(v00(),v01(),v02(),
        v10(),v11(),v12(),
        v20(),v21(),v22()));
  
    a_result.multiply(1./deter);  
  
    return true;
  }

private:
  T D33(const T& a_11,const T& a_12,const T& a_13,const T& a_21,const T& a_22,const T& a_23,const T& a_31,const T& a_32,const T& a_33) const {
    return (a_11 * (a_22 * a_33 - a_32 * a_23)- a_21 * (a_12 * a_33 - a_32 * a_13) + a_31 * (a_12 * a_23 - a_22 * a_13));
  }
};

template <class T>
inline sqm<T> operator-(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res(a1);
  res.subtract(a2);
  return res;
}
template <class T>
inline sqm<T> operator+(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res(a1);
  res.add(a2);
  return res;
}
template <class T>
inline sqm<T> operator*(const sqm<T>& a1,const sqm<T>& a2) {
  sqm<T> res;
  a1.mx_mul(a2,res);
  return res;
}

template <class T>
inline sqm<T> commutator(const sqm<T>& a1,const sqm<T>& a2) {
  return a1*a2-a2*a1;
}
template <class T>
inline sqm<T> anticommutator(const sqm<T>& a1,const sqm<T>& a2) {
  return a1*a2+a2*a1;
}

template <class T>
class cbm : public array<T> { //cbm for cubic matrix.
public:
  cbm():array<T>(3,4){}
  virtual ~cbm() {}
public:
  cbm(const cbm& a_from):array<T>(a_from){}
  cbm& operator=(const cbm& a_from){
    inlib::array<T>::operator=(a_from);
    return *this;
  }
public:
  T value(unsigned int a1,unsigned int a2,unsigned int a3) const { 
    return array<T>::m_vector[a1 + a2 * 4 + a3 * 16];
  }

  void v(unsigned int a1,unsigned int a2,unsigned int a3,const T& a_value) { 
    array<T>::m_vector[a1 + a2 * 4 + a3 * 16] = a_value;
  }

  void v000(const T& a_value){array<T>::m_vector[0+0*4+0*16] = a_value;}
  void v100(const T& a_value){array<T>::m_vector[1+0*4+0*16] = a_value;}
  void v200(const T& a_value){array<T>::m_vector[2+0*4+0*16] = a_value;}
  void v300(const T& a_value){array<T>::m_vector[3+0*4+0*16] = a_value;}
  void v010(const T& a_value){array<T>::m_vector[0+1*4+0*16] = a_value;}
  void v110(const T& a_value){array<T>::m_vector[1+1*4+0*16] = a_value;}
  void v210(const T& a_value){array<T>::m_vector[2+1*4+0*16] = a_value;}
  void v310(const T& a_value){array<T>::m_vector[3+1*4+0*16] = a_value;}
  void v020(const T& a_value){array<T>::m_vector[0+2*4+0*16] = a_value;}
  void v120(const T& a_value){array<T>::m_vector[1+2*4+0*16] = a_value;}
  void v220(const T& a_value){array<T>::m_vector[2+2*4+0*16] = a_value;}
  void v320(const T& a_value){array<T>::m_vector[3+2*4+0*16] = a_value;}
  void v030(const T& a_value){array<T>::m_vector[0+3*4+0*16] = a_value;}
  void v130(const T& a_value){array<T>::m_vector[1+3*4+0*16] = a_value;}
  void v230(const T& a_value){array<T>::m_vector[2+3*4+0*16] = a_value;}
  void v330(const T& a_value){array<T>::m_vector[3+3*4+0*16] = a_value;}

  void v001(const T& a_value){array<T>::m_vector[0+0*4+1*16] = a_value;}
  void v101(const T& a_value){array<T>::m_vector[1+0*4+1*16] = a_value;}
  void v201(const T& a_value){array<T>::m_vector[2+0*4+1*16] = a_value;}
  void v301(const T& a_value){array<T>::m_vector[3+0*4+1*16] = a_value;}
  void v011(const T& a_value){array<T>::m_vector[0+1*4+1*16] = a_value;}
  void v111(const T& a_value){array<T>::m_vector[1+1*4+1*16] = a_value;}
  void v211(const T& a_value){array<T>::m_vector[2+1*4+1*16] = a_value;}
  void v311(const T& a_value){array<T>::m_vector[3+1*4+1*16] = a_value;}
  void v021(const T& a_value){array<T>::m_vector[0+2*4+1*16] = a_value;}
  void v121(const T& a_value){array<T>::m_vector[1+2*4+1*16] = a_value;}
  void v221(const T& a_value){array<T>::m_vector[2+2*4+1*16] = a_value;}
  void v321(const T& a_value){array<T>::m_vector[3+2*4+1*16] = a_value;}
  void v031(const T& a_value){array<T>::m_vector[0+3*4+1*16] = a_value;}
  void v131(const T& a_value){array<T>::m_vector[1+3*4+1*16] = a_value;}
  void v231(const T& a_value){array<T>::m_vector[2+3*4+1*16] = a_value;}
  void v331(const T& a_value){array<T>::m_vector[3+3*4+1*16] = a_value;}

  void v002(const T& a_value){array<T>::m_vector[0+0*4+2*16] = a_value;}
  void v102(const T& a_value){array<T>::m_vector[1+0*4+2*16] = a_value;}
  void v202(const T& a_value){array<T>::m_vector[2+0*4+2*16] = a_value;}
  void v302(const T& a_value){array<T>::m_vector[3+0*4+2*16] = a_value;}
  void v012(const T& a_value){array<T>::m_vector[0+1*4+2*16] = a_value;}
  void v112(const T& a_value){array<T>::m_vector[1+1*4+2*16] = a_value;}
  void v212(const T& a_value){array<T>::m_vector[2+1*4+2*16] = a_value;}
  void v312(const T& a_value){array<T>::m_vector[3+1*4+2*16] = a_value;}
  void v022(const T& a_value){array<T>::m_vector[0+2*4+2*16] = a_value;}
  void v122(const T& a_value){array<T>::m_vector[1+2*4+2*16] = a_value;}
  void v222(const T& a_value){array<T>::m_vector[2+2*4+2*16] = a_value;}
  void v322(const T& a_value){array<T>::m_vector[3+2*4+2*16] = a_value;}
  void v032(const T& a_value){array<T>::m_vector[0+3*4+2*16] = a_value;}
  void v132(const T& a_value){array<T>::m_vector[1+3*4+2*16] = a_value;}
  void v232(const T& a_value){array<T>::m_vector[2+3*4+2*16] = a_value;}
  void v332(const T& a_value){array<T>::m_vector[3+3*4+2*16] = a_value;}

  void v003(const T& a_value){array<T>::m_vector[0+0*4+3*16] = a_value;}
  void v103(const T& a_value){array<T>::m_vector[1+0*4+3*16] = a_value;}
  void v203(const T& a_value){array<T>::m_vector[2+0*4+3*16] = a_value;}
  void v303(const T& a_value){array<T>::m_vector[3+0*4+3*16] = a_value;}
  void v013(const T& a_value){array<T>::m_vector[0+1*4+3*16] = a_value;}
  void v113(const T& a_value){array<T>::m_vector[1+1*4+3*16] = a_value;}
  void v213(const T& a_value){array<T>::m_vector[2+1*4+3*16] = a_value;}
  void v313(const T& a_value){array<T>::m_vector[3+1*4+3*16] = a_value;}
  void v023(const T& a_value){array<T>::m_vector[0+2*4+3*16] = a_value;}
  void v123(const T& a_value){array<T>::m_vector[1+2*4+3*16] = a_value;}
  void v223(const T& a_value){array<T>::m_vector[2+2*4+3*16] = a_value;}
  void v323(const T& a_value){array<T>::m_vector[3+2*4+3*16] = a_value;}
  void v033(const T& a_value){array<T>::m_vector[0+3*4+3*16] = a_value;}
  void v133(const T& a_value){array<T>::m_vector[1+3*4+3*16] = a_value;}
  void v233(const T& a_value){array<T>::m_vector[2+3*4+3*16] = a_value;}
  void v333(const T& a_value){array<T>::m_vector[3+3*4+3*16] = a_value;}

};

}}

#endif