vsr_generic_op.h

Go to the documentation of this file.

#ifndef VSR_GENERIC_OP_H_INCLUDED
#define VSR_GENERIC_OP_H_INCLUDED   

#include "vsr_multivector.h"
#include "../util/vsr_constants.h" 
#include <vector>  

namespace vsr {      
 
 

/*-----------------------------------------------------------------------------
 *  Group Definitions for Doxygen...
 *-----------------------------------------------------------------------------*/
  
namespace nga{ 
    
template<int DIM>
struct Proj{
  
  typedef NEVec<DIM> TVec;       
  typedef NEVec<DIM-1> OneDown;  
     
  static auto Call( VSR_PRECISION dist, const TVec& v ) RETURNS (
    ( Proj<DIM-1>::Call( dist, v.template cast<OneDown>() * ( dist / (dist - v[DIM-1] ) ) ) )
  )  
    
  template<int DIM2>
  static auto Ortho( const TVec& v ) RETURNS (
    ( v.template cast<NEVec<DIM2> >() )
  )   

  static auto Ortho3( const TVec& v ) RETURNS (
    ( v.template cast<NEVec<3>>() )
  )  
  
 
  
  static VSR_PRECISION Val( VSR_PRECISION dist, const TVec & v ) { 
    return dist / (dist - v[DIM-1] )  * Proj<DIM-1>::Val(dist, OneDown(v) ); 
  }
  // static VSR_PRECISION Val( VSR_PRECISION dist, const Vec& v) {
  //   return Proj< DIM, 4, DIM==TARGET >::Call(dist, v);
  // } 
};    
     
template<>
struct Proj<3>{  
  typedef NEVec<3> TVec;
  static TVec Call(VSR_PRECISION dist, const TVec& v) { return v; }  
  static VSR_PRECISION Val(VSR_PRECISION dist, const TVec & v ) { return 1.0; }
   static TVec Ortho3( const TVec& v ) { return v; }
}; 

//lift up to 3 from 2
template<>
struct Proj<2>{  
  typedef NEVec<3> TVec;
  static TVec Call(VSR_PRECISION dist, const TVec& v) { return v; }  
  static VSR_PRECISION Val(VSR_PRECISION dist, const TVec & v ) { return 1.0; }
   static TVec Ortho3( const TVec& v ) { return v; }
}; 
 

 template< class X>
 constexpr VSR_PRECISION dot(X x){
  return x*x;
}

 template<class X, class ... XS>
 constexpr VSR_PRECISION dot(X x, XS...xs){
  return (x*x) + dot(xs...);
 }  

 
struct Op{  
  
  template<class A>  
  static auto dual( const A& a ) RETURNS (
    a.dual()
  ) 

  template<class A>  
  static auto undual( const A& a ) RETURNS (
    a.undual()
  )

  template<class A>  
  static auto duale( const A& a ) RETURNS (
    a.duale()
  )
   
  template<class A>  
  static auto unduale( const A& a ) RETURNS (
    a.unduale()
  )
  
  template<class T> static auto dl( const T& t ) RETURNS ( dual(t) )
  template<class T> static auto udl( const T& t) RETURNS ( udual(t) )
  template<class T> static auto dle( const T& t) RETURNS (duale(t))
  template<class T> static auto udle( const T& t) RETURNS (unduale(t))
 
  
  template<class A>
  static constexpr bool sign(const A& a, const A& b) {
    return (a / b)[0] > 0 ? 1 : 0;
  } 

  template<class A> static constexpr bool sn(const A& a, const A& b){
    return sign(a,b);
  }

  template< class A, class B>
  static constexpr auto project(const A& a, const B& b) RETURNS ( 
    (a <= b ) / b
  ) 

  
  template< class A, class B>
  static constexpr auto reject(const A& a, const B& b) RETURNS ( 
    (a ^ b ) / b
  )

  template< class A, class B>
  static constexpr auto pj(const A& a, const B& b) RETURNS ( project(a,b) )
  template< class A, class B>
  static constexpr auto rj(const A& a, const B& b) RETURNS ( reject(a,b) )

};



struct Euc{
    template<class algebra>
    GAVec< typename algebra::up > static hom(const GAVec<algebra>& v){
      using up = typename algebra::up;
      return GAVec<up>(v) + GAE<up>::e<up::Dim>(1.0);
    }

//    template< bits::type N, class T >
//    NEVec<N+1,T> hom(const NEVec<N,T>& v){
//      return NEVec<N+1,T>(v) + NE<T>::e<N+1>(1.0);
//    }


//    template<class algebra>
//    GAVec<typename algebra::> hom(const GAPnt< Algebra >& p){
//      return NEVec<N-1,T>(p[0],p[1],p[2]) + NE<T>::e<N-1>(1.0);
//    }
     
//    template<bits::type N, class T>
//    NEVec<N-1, T> hom(const NPnt<N, T>& p){
//      return NEVec<N-1,T>(p[0],p[1],p[2]) + NE<T>::e<N-1>(1.0);
//    }
};



 struct Gen{


    template<class A>
    static auto rot ( const A& b ) -> decltype( b + 1 ) {  
      //  printf("me");
      VSR_PRECISION  c = sqrt(- ( b.wt() ) );
      VSR_PRECISION sc = -sin(c);
      if (c != 0) sc /= c;
      return b * sc + cos(c);
    }
    
    template<class A> 
    static auto rotor( const A& b ) -> decltype( b + 1 ) {
      return rot(b); 
    }

         
     template<class algebra>
     static auto log(const GARot<algebra>& r) -> GABiv<algebra> {
                 
        using TBiv = GABiv< algebra >;

        VSR_PRECISION t = r.template get<0>();                           //<--- Field Value from Rotor, or r[0]

        TBiv b = r.template cast<TBiv>();

        VSR_PRECISION n = b.rnorm();

        if (n <= 0) {
            if (t < 0) {
                return TBiv(PI);
            } else {
                return TBiv(); 
            }
        }

        VSR_PRECISION s = atan2( n, t );
        return b * ( s / n);
    } 
  
  
     template<class algebra>
     static auto log(const GABst<algebra>& r) -> GAPar<algebra> {
         
        using TPar = GAPar<algebra>;  
          
        VSR_PRECISION n;

        TPar p;
        p = r; //extract 2-blade part
        VSR_PRECISION td = p.wt(); //get scalar

        if (td > 0 ) { VSR_PRECISION s2 = sqrt(td);  n = asinh( s2 ) / s2; }
        else if ( td == 0 ) { n = 1; }
        else if (td < 0 ) { VSR_PRECISION s2 = sqrt(-td); n = atan2(s2, r[0] ) / s2; }

        return p * n;
  
      }
  
    template< class A >
    static auto pl( const GARot<A>& r) -> GABiv<A> { 
          using TBiv =  GABiv<A>;
          TBiv b = r.template cast<TBiv>();
          VSR_PRECISION t = b.rnorm(); // use rnorm or norm here?
          if (t == 0 ) return TBiv(1);
          return b / t;
      } 

     template<class A> static auto plane( const A& r) -> GABiv<A>
     { return pl(r); }

     template< class A > 
     static VSR_PRECISION iphi( const GARot<A>& r) { 
          using TBiv = GABiv<A>;  
          return TBiv ( log(r) * -2 ).norm();
      }

     template< class A > 
     static auto aa (const GARot<A>& r) -> GARot<A> {
         using TRot = GARot<A>;
         using TVec = GAVec<A>;
      
          TVec v = Op::dle( pl( r ) ) ;    
          VSR_PRECISION deg = iphi(r) * ( -180 / PI );

          return TRot(deg, v[0], v[1], v[2]);
      }
      
      template< class A >
      static auto axisAngle( const A& r) RETURNS ( aa(r) )
  
     
     template<class A, class B>
     static auto trv ( const Multivector<A,B>& a ) RETURNS (
       a.template copy< GATnv<A> >() + 1
     )

    template< class A > static auto transversor( const A& a ) RETURNS ( trv(a) ) 
    
    
      /*
          Generate transversor  as exponential of a direction vector
          @param a tangent vector
      */

  //  template< typename ... T >
  //  auto trv ( T ... v ) ->  NTrv<sizeof...(T)+2> {      
  //    return ( NTnv<sizeof...(T)+2>(v...) ) + 1;
  //  }
   
    template< typename ... T> static auto transversor( T ... v ) RETURNS ( trv(v...) )
    
    template<class A, class B>
    static GATrs<A> trs ( const Multivector<A,B>& a )  {
        return (a.template copy< GADrv<A> > () * -.5 ) + 1;
    }
    template<class A> static auto translator (const A& a) RETURNS ( trs(a) ) 
    
    
    template<  typename ... Ts >
    static auto trs ( Ts ... v ) ->  NTrs<sizeof...(Ts)+2> {      
      return ( NDrv<sizeof...(Ts)+2>(v...) * -.5 ) + 1;
    } 
 //   template<typename ... T> auto translator(T ... v) RETURNS ( trs( v...) )
    
      template<bits::type N, class T>
      static constexpr auto dil( T t) -> NDil<N> {
          return NDil<N>( cosh( t *.5 ), sinh( t * .5 ) );
      } 
      template<class T> static constexpr auto dilator( T t ) RETURNS ( dil(t) )

      
      template<class A, class T>
      static auto dil(const GAPnt<A>& p, T t) -> GATsd<A>  {
          return GATsd<A>( GADil<A>( cosh( t*.5 ), sinh( t*.5 ) ) ).trs( p );
      }
      template< class A, class T > static constexpr auto dilator ( const A& p, T t) RETURNS ( dil(p,t) ) 
    
     
      template<class A>  
      static auto bst(const GAPar<A>& tp) -> decltype( tp + 1 ) { 

          VSR_PRECISION norm; VSR_PRECISION sn; VSR_PRECISION cn;

          VSR_PRECISION td = tp.wt(); 

          if (td < 0) { norm =  sqrt( - td );  sn = -sin(norm) / norm; cn = cos(norm); } //note, changed to cos from cosh
          else if (td > 0) { norm = sqrt(td); sn = -sinh(norm) / norm; cn = cosh(norm); }
          else if (td == 0) { norm = 0; sn = -1; cn = 1; }

          return (tp * sn) + cn;
    }
    template<class A> static auto boost( const A& a ) RETURNS ( bst(a) )
    
    // /*!
    //           Generate general rotation as exponential of anonymous 
    //           Implemented from "Square Root and Logarithm of Rotors. . ." by Dorst and Valkenburg, 2011
    //           @param Point Pair generator
    //       */ 
    // template<bits::type DIM, class A>  
    //       auto gen(const CGAMultivector<DIM,A>& tp) -> decltype( tp + 1 ) { 
    // 
    //   VSR_PRECISION norm; VSR_PRECISION sn; VSR_PRECISION cn;
    // 
    //           VSR_PRECISION td = tp.wt(); 
    // 
    //           if (td < 0) { norm =  sqrt( - td );  sn = -sin(norm) / norm; cn = cosh(norm); }
    //           else if (td > 0) { norm = sqrt(td); sn = -sinh(norm) / norm; cn = cosh(norm); }
    //           else if (td == 0) { norm = 0; sn = -1; cn = 1; }
    // 
    //           return (tp * sn) + cn;
    // } 
    

    
     //feed in vectors!  
     template<class A> 
     static auto ratio( 
         const GAVec<A>& a, 
         const GAVec<A>& b ) -> decltype( (a*b) ) {
          
         using TVec = GAVec<A>;//typename NVec<DIM>::Space::Vec;
         using TBiv = GABiv<A>;//typename NVec<DIM>::Space::Biv;
         using TRot = GARot<A>;//decltype( (a^b) + 1);
         
         VSR_PRECISION s = ( a <= b )[0];              
         //180 degree check
         if ( a == b.conjugation() ){//fabs ((a<=b)[0]) > .999999) {//a == b.conjugation() ) {
           //printf("180!\n");
           if ( a == TVec::y || a == -TVec::y ) {
             return rot( TBiv::xy * PIOVERTWO );
           }        
           return rot( a ^ TVec::y * PIOVERTWO); //mind the ordering of blades
         }         
      
         VSR_PRECISION ss = 2 * (s+1);
         //VSR_PRECISION n = ( ss >= 0 ? sqrt ( ss ) : 0 );
         VSR_PRECISION n = ( ss >= 0 ? sqrt ( ss ) : -sqrt(-ss) );

         TRot r = ( b * a ) ; //cout << r << endl;
         r[0] += 1;  
         if (!FERROR(n)) r /= n;
         if ( r == TRot() ) {  return TRot(1); }//else cout << r << endl; //printf("0 in Gen::ratio\n");
         return r;    
      } 

 }; 
   



 struct Round {                          

  template< class A, class B >
  static constexpr GAPnt<A> 
  null( const Multivector<A,B>& v ) {  
       using TVec = GAVec<A>;
       using TOri = GAOri<A>;
       using TInf = GAInf<A>;
       return v.template copy< TVec >() + TOri( 1 ) + TInf( v.template copy< TVec >().wt() / 2.0 );
   } 
    
 template< class ... T>
 static constexpr NPnt< sizeof...(T) + 2>
 null( T ... v) { 
   using TVEC = NVec<sizeof...(T)+2>;
   using TORI = NOri<sizeof...(T)+2>;
   using TINF = NInf<sizeof...(T)+2>;
   
   return TVEC(v...) + TORI(1) + TINF( TVEC(v...).wt() / 2.0 ); 
 }

  template< class ... T>
  static constexpr NPnt< sizeof...(T) + 2>
  point( T ... v) { 
    using TVEC = NVec<sizeof...(T)+2>;
    return null( TVEC(v...) );
  }


  template< class ... T >
  static auto 
  dls(VSR_PRECISION r, T ... v ) ->  NPnt< sizeof...(T) + 2 >  { 
    using TPNT = NPnt< sizeof...(T) + 2 >;
        TPNT s = point(v...);
        ( r > 0) ? s.template get< bits::infinity< sizeof...(T) + 2>() >() -= .5 * (r * r) 
                 : s.template get< bits::infinity< sizeof...(T) + 2 >() >() += .5 * (r*r);
    return s;     
    } 
    template< class ... T > static auto dualSphere( VSR_PRECISION r, T ... v ) RETURNS ( dls(r, v...) )

    
    template<class A, class B >
    static auto 
    dls( const Multivector<A,B>& v, VSR_PRECISION r = 1.0 ) -> GAPnt<A>{
        
        auto s = null(v);
        ( r > 0) ? s.template get< bits::infinity< A::dim >() >() -= .5 * (r * r)
                 : s.template get< bits::infinity< A::dim >() >() += .5 * (r*r);
        return s;     
    }
    
    template< class T > static auto dualSphere( const T& t, VSR_PRECISION r = 1.0 ) RETURNS ( dls(t, r) )


    template< class S >
    static auto 
    sphere( const S& v, VSR_PRECISION r = 1.0 ) RETURNS(
      dls( v, r )    
    )



    template< class A >
    static GADls<A> 
    dls_pnt( const GAPnt<A>& p, VSR_PRECISION r = 1.0 ) {
        GAPnt<A> s = p;
        (r > 0) ? s.template get< bits::infinity<A::dim>() >() -= .5 * (r * r) 
                : s.template get< bits::infinity<A::dim>() >() += .5 * (r*r);
        return s;
    } 

    template<class A, class B>
    static constexpr GAPnt<A>
    center( const Multivector<A,B>& s) {
        return  ( s  / ( GAInf<A>(-1) <= s ) ).template cast< GAPnt<A> >();
    }

    template<class A, class B>
    static constexpr GAPnt<A>
    cen( const Multivector<A,B>& s) {
        return center(s);
    }


    template<class A>
    static constexpr typename A::space::point 
    location(const A& s) {
      return null ( cen ( s ) ); 
    }
    
    template<class A>
    static constexpr typename A::space::point
    loc(const A& s) { return location(s); }   
  

    template<class A> 
    static VSR_PRECISION 
    size( const A& r, bool dual){
        auto s = typename A::space::infinity(1) <= r;
        return ( ( r * r.inv() ) / ( s * s ) * ( (dual) ? -1.0 : 1.0 )  )[0];
    } 
  template<class T> 
    static constexpr VSR_PRECISION 
    radius( const T& s ){
        return sqrt ( fabs ( size(s, false) ) );
    } 
    template<class T> static constexpr VSR_PRECISION rad( const T& t) { return radius(t); }
   
    template<class A>
    static VSR_PRECISION 
    curvature(const A& s){
        VSR_PRECISION r = rad( s );     
        return (r==0) ? 10000 : 1.0 / rad(s);
    }
    
    template<class T> 
    static constexpr VSR_PRECISION 
    cur( const T& t) { return curvature(t); }


    template<class A>
    static constexpr VSR_PRECISION 
    dsize( const GAPnt<A>& dls ){
        return (dls * dls)[0];
    }
   
  template<class A>
    static constexpr VSR_PRECISION 
    squaredDistance(const GAPnt<A>& a, const GAPnt<A> b){
        return ( (a <= b)[0] ) * -2.0;
    }
    template< class A> static constexpr VSR_PRECISION sqd( const A& a, const A& b) { return squaredDistance(a,b); }

    
  template<class A>
    static constexpr VSR_PRECISION 
    distance(const GAPnt<A>& a, const GAPnt<A> b){
        return sqrt( fabs(sqd(a,b) ) );
    }
    template< class A> static constexpr VSR_PRECISION dist( const A& a, const A& b) { return distance(a,b); }

    template<class A>
    static std::vector< GAPnt<A> > 
    split(const GAPar<A>& pp){
        
      std::vector< GAPnt<A> > pair;
          
      VSR_PRECISION r = sqrt( fabs( ( pp <= pp )[0] ) );
          
      //dual line in 2d, dual plane in 3d
      auto d = GAInf<A>(-1) <= pp;

      GABst<A> bstA; bstA = pp;        
      GABst<A> bstB; bstB = pp;        
          
      bstA += GASca<A>(r);
      bstB -= GASca<A>(r);
                  
      GAPnt<A> pA;
      pA = ( bstA ) / d;
      GAPnt<A> pB;
      pB = ( bstB ) / d;
                  
      pair.push_back(pA);
      pair.push_back(pB);

      return pair;
    }

    template<class A>
    static std::vector< GAPnt<A> > 
    splitLocation(const GAPar<A>& pp){  
      auto tp = split(pp);
      for (auto& i : tp) i = location(i);
      return tp;
    }
   
    template<class A> 
    static GAPnt<A> 
    split(const GAPar<A>& pp, bool bFirst){
        
        VSR_PRECISION r = sqrt( fabs( ( pp <= pp )[0] ) );
        
        auto d = GAInf<A>(-1) <= pp;
        
        GABst<A> bst = pp + ( bFirst ? r : -r ); 
        
        return ( ( bst ) / d ).template cast< GAPnt<A> >();  

    }  

    template<class A>
    static std::vector< GAPnt<A> >
    split( const GACir<A>& nc ){
      return split( nc.dual() );
    }

    template<class A>
    static constexpr auto 
    direction( const A& s ) RETURNS(
        ( ( typename A::space::infinity(-1) <= s ) ^ typename A::space::infinity(1) )
    ) 
    
    template<class A> static constexpr auto dir( const A&s ) RETURNS( direction(s) ) 


    template<class A>
    static constexpr auto 
    carrier(const A& s) RETURNS(
        s ^ typename A::space::infinity(1)
    ) 
    template<class A> 
    static constexpr auto car( const A&s ) RETURNS( carrier(s) ) 
    
    template<class A>
    static constexpr typename A::space::dual_sphere
    surround( const A& s) {
        return typename A::space::dual_sphere( s / ( s ^ typename A::space::infinity(1) ));
    } 
    
    template<class A> 
    static constexpr auto 
    sur( const A& s) RETURNS ( surround(s) )

    
     template<class A, class S>
     static constexpr auto 
     produce(const A& dls, const S& flat) RETURNS (
         dls ^ ( ( dls <= ( flat.inv() * typename A::space::infinity(1) ) )  * -1.0 ) 
     )

     template<class A>
     static constexpr auto
     real(const A& s) RETURNS (
         produce( 
                Round::dls( Round::loc( s ), -Round::rad( Round::sur( s ) ) ), 
                typename A::space::origin(-1) <= Round::dir( s ) 
              )
     )


     template<class A>
     static constexpr auto
     imag(const A& s) RETURNS (
         produce( 
                Round::dls( Round::loc( s ), Round::rad( Round::sur( s ) ) ), 
                typename A::space::origin(-1) <= Round::dir( s ) 
              )
     )
     template<class A>
     static constexpr GADls<A> 
      at( const GADls<A> & c, const GADls<A>& p) {
        return GADls<A>( p <= (c^GAInf<A>(1) ) );
     }
    
//  /*!
//  Direct Round From coordinates and Euclidean Bivector
//  Note: round will be imaginary if dual sphere is real . . .
//  */  
//  template<bits::type DIM, class B, class A>
//  auto round(const CGAMultivector<DIM, B>& s, const A& flat, VSR_PRECISION radius=1.0) RETURNS (
//     round( dls(s, radius*-1.0), flat ); 
//  ) 
  
 
  

    template<class A>
     static constexpr GAPnt<A> pnt(const GADls<A> & dls, const GAVec<A>& flat){
        return split( produce(dls, flat), true ); // cout << "y" << endl; 
     }
   // template<class A



    
    template<class A>
    static GAVec<A> 
    vec(const GACir<A>& c, VSR_PRECISION theta = 0){ 
        using TBIV = GABiv<A>;
    
        GADll<A> axis = (GAInf<A>(1) <= c).runit();        
        return GAVec<A>::x.sp( nga::Gen::ratio( TBIV::xz.duale(), TBIV(axis).duale() ) * nga::Gen::rot(TBIV::xz * theta) );    
    }    


    template<class A>
    static GAPar<A> 
    par_cir(const GACir<A>& c, VSR_PRECISION t) {  
        using TBIV = GABiv<A>;
        using TVEC = GAVec<A>;
    
        //NDll<DIM> axis = (NInf<DIM>(-1) <= c).runit();      
        //NRot<DIM> rot = Gen::ratio( TBIV::xz.duale(), TBIV(axis).duale() );
       // NVec<DIM> vec = NVec<DIM>::x.sp( rot * Gen::rot(TBIV::xz * t/2.0)); //BIG CHANGE . . .
        auto normal = TVEC( Round::carrier(c).dual() ).unit();
        auto rot = nga::Gen::ratio( TVEC::z, normal );
        auto vec = TVEC::x.spin( nga::Gen::rot( TBIV::xy.spin( rot ) * t/2.0 ) * rot ); 
                    
        return produce( sur( c ), vec );
   }       
  
    template<class A>
    static GAPnt<A> 
    pnt_cir(const GACir<A>& c, VSR_PRECISION t) { 
      return null( split( par_cir(c,t), true) );
    }


    /* template<bits::type DIM> */
    /* NPnt<DIM> */ 
    /* pointOnCircle(const NCir<DIM>& c, VSR_PRECISION t){ */
    /*   return null( split( par_cir(c,t), true) ); */
   /* } */



  /* template<bits::type DIM> */
  /* NPnt<DIM> pnt_dls( */


 };
       

  struct Flat {          

      template<class A,class B> 
      static constexpr auto direction( const Multivector<A,B>& f) RETURNS(
          GAInf<A>(-1) <= f
      )  
    
      template<class A,class B> 
      static constexpr auto dir( const Multivector<A,B>& f) RETURNS(
          direction(f)
      )  
    
    
      template<class A>
      static constexpr typename A::space::Pnt location(const A& f, const typename A::space::Pnt& p, bool dual){
          using TPnt = typename A::space::Pnt; 
          return dual ? TPnt( ( p ^ f ) / f ) :  TPnt ( ( p <= f ) / f );
      } 

      
      template<class A, class P>
      static constexpr P loc(const A& f, const P& p, bool dual){
          return location(f,p,dual);
      }

     template<class A>
     static constexpr typename A::value_t wt(const A& f, bool bDual){
         using TOri = typename A::space::origin;
         return bDual ? ( TOri(1) <= dir( f.undual() ) ).wt() : ( TOri(1) <= dir(f) ).wt();
     } 
     template<class A>
     static constexpr auto dlp( const GAPnt<A>& pnt, const GADrv<A>& drv) RETURNS (
        pnt <= drv
    )      

      
      template<typename ...T>
      static constexpr NLin<sizeof...(T)+2> line( T ... v ){
        return nga::Round::null( NVec<sizeof...(T)+2>() ) ^ nga::Round::null(v...) ^ NInf<sizeof...(T)+2>(1);
      }

      template<typename ...T>
      static constexpr NCir<sizeof...(T)+2> dline( T ... v ){
        return nga::Round::null( NVec<sizeof...(T)+2>() ) ^ nga::Round::null(v...) ^ nga::Round::dls(NVec<sizeof...(T)+2>(), 1.0 );
      }


  };
   


   struct Tangent {

//    /*! Direction of Tangent Element (similar formulation to Rounds) 
//        @ingroup direction
//
//        @param a Direct Tangent Element
//        @return @ref direction type
//    */
//    template <class A>
//    static constexpr auto dir (const A& a) RETURNS (
//        ( typename A::space::infinity(-1) <= a ) ^ typename A::space::infinity(1)
//    )
//
    template< class A >
    static constexpr typename A::space::point loc( const A& s){
        return ( s / typename A::space::infinity(-1) <= s );
    }

    template< class A >
    static constexpr auto at( const A& r, const typename A::space::point& p) RETURNS(
       p <= r.inv()
    )

    template<class A>
    static typename A::value_t wt(const A& s){
        using TOri = typename A::space::origin;
        return ( TOri(1) <= Round::dir(s) ).wt();
    }
      
   };

}//nga::
  
  //------------------------------------------
  
 //METHODS (MOTORS IMPLEMENTED SEPARATELY, IN SPECIFIC INSTANTIATIONS)
 template<class Algebra, class B>
 Multivector<Algebra, typename Algebra::vector_basis> Multivector<Algebra,B>::null() const{
   return nga::Round::null(*this);
 } 
  
 template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::rot(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::rot(t) );
 } 
 template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::rotate(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::rot(t) );
 } 

 template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::trs(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::trs(t) );
 } 
 template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::translate(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::trs(t) );
 } 


 template<class Algebra, class B> template<class ... Ts>
 Multivector<Algebra,B> Multivector<Algebra,B>::trs(Ts ... v) const{
   return this->sp( nga::Gen::trs( GADrv<Algebra>(v...) ) );
 } 
  template<class Algebra, class B> template<class ... Ts>
 Multivector<Algebra,B> Multivector<Algebra,B>::translate(Ts ... v) const{
   return this->sp( nga::Gen::trs( GADrv<Algebra>(v...) ) );
 } 

  template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::trv(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::trv(t) );
 } 
   template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::transverse(const Multivector<Algebra,A>& t) const{
   return this->sp( nga::Gen::trv(t) );
 } 

 template<class Algebra, class B> template<class ... Ts>
 Multivector<Algebra,B> Multivector<Algebra,B>::trv(Ts ... v) const{
   return this->sp( nga::Gen::trv( GATrv<Algebra>(v...) ) );
 } 
  template<class Algebra, class B> template<class ... Ts>
 Multivector<Algebra,B> Multivector<Algebra,B>::transverse(Ts ... v) const{
   return this->sp( nga::Gen::trv( GATrv<Algebra>(v...) ) );
 } 
     
   template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::dil(const Multivector<Algebra,A>& s, VSR_PRECISION t) const{
   return this->sp( nga::Gen::dil(s,t) );
 } 
   template<class Algebra, class B> template<class A>
 Multivector<Algebra,B> Multivector<Algebra,B>::dilate(const Multivector<Algebra,A>& s, VSR_PRECISION t) const{
   return this->sp( nga::Gen::dil(s,t) );
 } 
 
 template<class Algebra, class B> template <class A>
Multivector<Algebra,B> Multivector<Algebra,B>::bst( const Multivector<Algebra,A>& t) const{
     return this -> sp ( nga::Gen::bst(t) );  
} 
template<class Algebra, class B> template <class A>
Multivector<Algebra,B> Multivector<Algebra,B>::boost( const Multivector<Algebra,A>& t) const{
     return this -> sp ( nga::Gen::bst(t) );  
} 



}  //vsr::

#endif