update vector2,3D

5 files changed, 797 insertions, 98 deletions

diff --git a/_test/meowpp.h b/_test/meowpp.h
index 86c44c1..8052f1a 100644
--- a/_test/meowpp.h
+++ b/_test/meowpp.h
@@ -1,6 +1,9 @@
 #ifndef   __meowpp_h__
 #define   __meowpp_h__
 
+#include "meowpp/geo/Vector2D.h"
+#include "meowpp/geo/Vector3D.h"
+
 #include "meowpp/Usage.h"
 #include "meowpp/oo/Register_Implement.h"
 
diff --git a/meowpp/geo/Vector2D.h b/meowpp/geo/Vector2D.h
index 7de3523..d112f57 100644
--- a/meowpp/geo/Vector2D.h
+++ b/meowpp/geo/Vector2D.h
@@ -1,66 +1,73 @@
 #ifndef   geo_Vector2D_H__
 #define   geo_Vector2D_H__
 
+#include "../math/Matrix.h"
+
 #include <cmath>
-#include "../math/utility.h"
 
 namespace meow{
   template<class Scalar>
   class Vector2D{
     private:
-      Scalar _x, _y, _w;
+      Scalar _x, _y;
     public:
+      Vector2D();
+      Vector2D(Vector2D const& __v);
+      Vector2D(Scalar const& __x, Scalar const& __y);
+      ~Vector2D();
+      //
+      Vector2D& copy(Vector2D const& __v);
+      //
+      Scalar const& x() const;
+      Scalar const& y() const;
+      Scalar const& x(Scalar const& __s);
+      Scalar const& y(Scalar const& __s);
+      Vector2D&    xy(Scalar const& __x, Scalar const& __y);
+      //
+      Vector2D positive() const;
+      Vector2D negative() const;
+      Vector2D right   () const;
       //
-      Vector2D(): _x(0), _y(0), _w(0){ }
-      Vector2D(Scalar const& x,
-               Scalar const& y)  : _x(   x), _y(   y), _w(   0){ }
-      Vector2D(Scalar const& x,
-               Scalar const& y,
-               Scalar const& w)  : _x(   x), _y(   y), _w(   w){ }
-      Vector2D(Vector2D const& v): _x(v._x), _y(v._y), _w(v._w){ }
+      Vector2D add(Vector2D const& __v) const;
+      Vector2D sub(Vector2D const& __v) const;
+      Vector2D mul(Scalar   const& __s) const;
+      Vector2D div(Scalar   const& __s) const;
       //
-      Scalar const& x() const{ return _x; }
-      Scalar const& y() const{ return _y; }
-      Scalar const& w() const{ return _w; }
-      Scalar      & x(Scalar const& s){ _x = s; return _x; }
-      Scalar      & y(Scalar const& s){ _y = s; return _y; }
-      Scalar      & w(Scalar const& s){ _w = s; return _w; }
-      Vector2D& operator()(Scalar const& x, Scalar const& y)                 { _x = x; _y = y;         return *this; }
-      Vector2D& operator()(Scalar const& x, Scalar const& y, Scalar const& w){ _x = x; _y = y; _w = w; return *this; }
-      Scalar&       operator[](size_t n)      { return (n == 0 ? _x : (n == 1 ? _y : _w)); }
-      Scalar const& operator[](size_t n) const{ return (n == 0 ? _x : (n == 1 ? _y : _w)); }
+      Scalar dot  (Vector2D const& __v) const;
+      Scalar cross(Vector2D const& __v) const;
+      Scalar length () const;
+      Scalar length2() const;
       //
-      Vector2D operator-() const{ return Vector2D(-_x, -_y, _w); }
-      Vector2D operator~() const{ return Vector2D(-_y,  _x, _w); }
+      Vector2D& normalize ();
+      Vector2D  normalized() const;
+      Vector2D& rotate    (Scalar const& __theta);
+      Vector2D  rotated   (Scalar const& __theta) const;
+      Vector2D& reflect   (Vector2D const& __v);
+      Vector2D  reflected (Vector2D const& __v) const;
       //
-      Vector2D operator+(Vector2D const& v) const{ return Vector2D(_x + v._x, _y + v._y, _w); }
-      Vector2D operator-(Vector2D const& v) const{ return Vector2D(_x - v._x, _y - v._y, _w); }
-      Vector2D operator*(Scalar   const& s) const{ return Vector2D(_x *    s, _y *    s, _w); }
-      Vector2D operator/(Scalar   const& s) const{ return Vector2D(_x /    s, _y /    s, _w); }
+      Matrix<Scalar> matrix(                    ) const;
+      Matrix<Scalar> matrix(Scalar const& __homo) const;
       //
-      Vector2D& operator+=(Vector2D const& v){ _x += v._x; _y += v._y; return *this; }
-      Vector2D& operator-=(Vector2D const& v){ _x -= v._x; _y -= v._y; return *this; }
-      Vector2D& operator*=(Scalar   const& s){ _x *=    s; _y *=    s; return *this; }
-      Vector2D& operator/=(Scalar   const& s){ _x /=    s; _y /=    s; return *this; }
+      Scalar const& operator()(size_t n) const;
+      Scalar&       operator()(size_t n);
+      Vector2D&     operator()(Scalar const& __x, Scalar const& __y);
       //
-      Scalar dot  (Vector2D const& v) const{ return _x * v._x + _y * v._y; }
-      Scalar cross(Vector2D const& v) const{ return _x * v._y - _y * v._x; }
-      Scalar length () const{ return sqrt(squ(_x) + squ(_y)); }
-      Scalar length2() const{ return      squ(_x) + squ(_y) ; }
+      Vector2D operator+() const;
+      Vector2D operator-() const;
+      Vector2D operator~() const;
       //
-      Vector2D normal() const{
-        Scalar len(length());
-        return Vector2D(_x / len, _y / len, _w);
-      }
-      Vector2D rotation(Scalar theta) const{
-        Vector2D i(cos(-theta), sin(-theta));
-        return Vector2D(i.dot(*this), i.cross(*this), _w);
-      }
-      Vector2D reflection(Vector2D const& v) const{
-        return v * v.dot(*this) * 2 / v.length2() - *this;
-      }
+      Vector2D operator+(Vector2D const& __v) const;
+      Vector2D operator-(Vector2D const& __v) const;
+      Vector2D operator*(Scalar   const& __s) const;
+      Vector2D operator/(Scalar   const& __s) const;
+      //
+      Vector2D& operator+=(Vector2D const& __v);
+      Vector2D& operator-=(Vector2D const& __v);
+      Vector2D& operator*=(Scalar   const& __s);
+      Vector2D& operator/=(Scalar   const& __s);
   };
-
 }
 
+#include "Vector2D.hpp"
+
 #endif // geo_Vector2D_H__
diff --git a/meowpp/geo/Vector2D.hpp b/meowpp/geo/Vector2D.hpp
new file mode 100644
index 0000000..e1d0f41
--- /dev/null
+++ b/meowpp/geo/Vector2D.hpp
@@ -0,0 +1,327 @@
+#include "Vector2D.h"
+
+#include "../math/utility.h"
+#include "../math/Matrix.h"
+
+#include <cmath>
+
+namespace meow{
+  template<class Scalar>
+  inline
+  Vector2D<Scalar>::Vector2D():
+  _x(0),
+  _y(0){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector2D<Scalar>::Vector2D(Vector2D const& __v):
+  _x(__v._x),
+  _y(__v._y){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector2D<Scalar>::Vector2D(Scalar const& __x, Scalar const& __y):
+  _x(__x),
+  _y(__y){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector2D<Scalar>::~Vector2D(){
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::copy(Vector2D const& __v){
+    _x = __v.x();
+    _y = __v.y();
+    return *this;
+  }
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector2D<Scalar>::x() const{
+    return _x;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector2D<Scalar>::y() const{
+    return _y;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector2D<Scalar>::x(Scalar const& __s){
+    _x = __s;
+    return x();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector2D<Scalar>::y(Scalar const& __s){
+    _y = __s;
+    return y();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::xy(Scalar const& __x, Scalar const& __y){
+    _x = __x;
+    _y = __y;
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::positive() const{
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::negative() const{
+    return Vector2D<Scalar>(-x(), -y());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::right() const{
+    return Vector2D<Scalar>(-y(), x());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::add(Vector2D const& __v) const{
+    return Vector2D<Scalar>(x() + __v.x(), y() + __v.y());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::sub(Vector2D const& __v) const{
+    return Vector2D<Scalar>(x() - __v.x(), y() - __v.y());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::mul(Scalar const& __s) const{
+    return Vector2D<Scalar>(x() * __s, y() * __s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::div(Scalar const& __s) const{
+    return Vector2D<Scalar>(x() / __s, y() / __s);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector2D<Scalar>::dot(Vector2D const& __v) const{
+    return (x() * __v.x() + y() * __v.y());
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector2D<Scalar>::cross(Vector2D const& __v) const{
+    return (x() * __v.y() - y() * __v.x());
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector2D<Scalar>::length() const{
+    return Scalar(sqrt(double(length2())));
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector2D<Scalar>::length2() const{
+    return squ(x()) + squ(y());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::normalize(){
+    return copy(normalized());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::normalized() const{
+    return div(length());
+  }
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::rotate(Scalar const& __theta){
+    return copy(rotated(__theta));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::rotated(Scalar const& __theta) const{
+    Scalar cs(cos(-double(__theta)));
+    Scalar sn(sin(-double(__theta)));
+    Vector2D<Scalar> new_x(cs, sn);
+    return Vector2D<Scalar>(new_x.dot(*this), new_x.cross(*this));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::reflect(Vector2D const& __v){
+    return copy(reflected(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::reflected(Vector2D const& __v) const{
+    return __v.mul(__v.dot(*this) * 2 / __v.length2()).sub(*this);
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Vector2D<Scalar>::matrix() const{
+    Matrix<Scalar> ret(2, 0, Scalar(0));
+    ret(0, 0) = x();
+    ret(1, 0) = y();
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Vector2D<Scalar>::matrix(Scalar const& __homo) const{
+    Matrix<Scalar> ret(3, 0, Scalar(0));
+    ret(0, 0) = x();
+    ret(1, 0) = y();
+    ret(2, 0) = __homo;
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector2D<Scalar>::operator()(size_t n) const{
+    return (n == 0 ? x() : y());
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar&
+  Vector2D<Scalar>::operator()(size_t n){
+    return (n == 0 ? _x : _y);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::operator()(Scalar const& __x, Scalar const& __y){
+    x(__x);
+    y(__y);
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator+() const{
+    return positive();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator-() const{
+    return negative();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator~() const{
+    return right();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator+(Vector2D const& __v) const{
+    return add(__v);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator-(Vector2D const& __v) const{
+    return sub(__v);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator*(Scalar const& __s) const{
+    return mul(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>
+  Vector2D<Scalar>::operator/(Scalar const& __s) const{
+    return div(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::operator+=(Vector2D const& __v){
+    return copy(add(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::operator-=(Vector2D const& __v){
+    return copy(sub(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::operator*=(Scalar const& __s){
+    return copy(mul(__s));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector2D<Scalar>&
+  Vector2D<Scalar>::operator/=(Scalar const& __s){
+    return copy(div(__s));
+  }
+}
+
diff --git a/meowpp/geo/Vector3D.h b/meowpp/geo/Vector3D.h
index 08b5705..6da1820 100644
--- a/meowpp/geo/Vector3D.h
+++ b/meowpp/geo/Vector3D.h
@@ -1,74 +1,77 @@
 #ifndef   geo_Vector3D_H__
 #define   geo_Vector3D_H__
 
+#include "../math/Matrix.h"
+
 #include <cmath>
-#include "../utility.h"
 
 namespace meow{
   template<class Scalar>
   class Vector3D{
     private:
-      Scalar _x, _y, _z, _w;
+      Scalar _x, _y, _z;
     public:
+      Vector3D();
+      Vector3D(Vector3D const& __v);
+      Vector3D(Scalar const& __x, Scalar const& __y, Scalar const& _z);
+      ~Vector3D();
+      //
+      Vector3D& copy(Vector3D const& __v);
+      //
+      Scalar const& x() const;
+      Scalar const& y() const;
+      Scalar const& z() const;
+      Scalar const& x(Scalar const& __s);
+      Scalar const& y(Scalar const& __s);
+      Scalar const& z(Scalar const& __s);
+      Vector3D&   xyz(Scalar const& __x, Scalar const& __y, Scalar const& __z);
       //
-      Vector3D(): _x(0), _y(0), _z(0), _w(0){ }
-      Vector3D(Scalar const& x,
-               Scalar const& y,
-               Scalar const& z)  : _x(   x), _y(   y), _z(   z), _w(   0){ }
-      Vector3D(Scalar const& x,
-               Scalar const& y,
-               Scalar const& z,
-               Scalar const& w)  : _x(   x), _y(   y), _z(   z), _w(   w){ }
-      Vector3D(Vector3D const& v): _x(v._x), _y(v._y), _z(v._z), _w(v._w){ }
+      Vector3D positive() const;
+      Vector3D negative() const;
+      Vector3D right   () const;
       //
-      Scalar const& x() const{ return _x; }
-      Scalar const& y() const{ return _y; }
-      Scalar const& z() const{ return _z; }
-      Scalar const& w() const{ return _w; }
-      Scalar      & x(Scalar const& s){ _x = s; return _x; }
-      Scalar      & y(Scalar const& s){ _y = s; return _y; }
-      Scalar      & z(Scalar const& s){ _z = s; return _z; }
-      Scalar      & w(Scalar const& s){ _w = s; return _w; }
-      Vector3D& operator()(Scalar const& x, Scalar const& y, Scalar const& z                 ){ _x = x; _y = y; _z = z;         return *this; }
-      Vector3D& operator()(Scalar const& x, Scalar const& y, Scalar const& z, Scalar const& w){ _x = x; _y = y; _z = z; _w = w; return *this; }
-      Scalar      & operator[](size_t n)      { return (n == 0 ? _x : (n == 1 ? _y : (n == 2 ? _z : _w))); }
-      Scalar const& operator[](size_t n) const{ return (n == 0 ? _x : (n == 1 ? _y : (n == 2 ? _z : _w))); }
+      Vector3D add(Vector3D const& __v) const;
+      Vector3D sub(Vector3D const& __v) const;
+      Vector3D mul(Scalar   const& __s) const;
+      Vector3D div(Scalar   const& __s) const;
       //
-      Vector3D operator-() const{ return Vector3D(-_x, -_y, -_z, _w); }
+      Scalar   dot  (Vector3D const& __v) const;
+      Vector3D cross(Vector3D const& __v) const;
+      Scalar length () const;
+      Scalar length2() const;
       //
-      Vector3D operator+(Vector3D const& v) const{ return Vector3D(_x + v._x, _y + v._y, _z + v._z, _w); }
-      Vector3D operator-(Vector3D const& v) const{ return Vector3D(_x - v._x, _y - v._y, _z - v._z, _w); }
-      Vector3D operator*(Scalar   const& s) const{ return Vector3D(_x *    s, _y *    s, _z *    s, _w); }
-      Vector3D operator/(Scalar   const& s) const{ return Vector3D(_x /    s, _y /    s, _z /    s, _w); }
+      Vector3D& normalize ();
+      Vector3D  normalized() const;
+      Vector3D& rotate    (Vector3D const& __axis, Scalar const& __theta);
+      Vector3D  rotated   (Vector3D const& __axis, Scalar const& __theta) const;
+      Vector3D& reflect   (Vector3D const& __v);
+      Vector3D  reflected (Vector3D const& __v) const;
       //
-      Vector3D& operator+=(Vector3D const& v){ _x += v._x; _y += v._y; _z += v._z; return *this; }
-      Vector3D& operator-=(Vector3D const& v){ _x -= v._x; _y -= v._y; _z -= v._z; return *this; }
-      Vector3D& operator*=(Scalar   const& s){ _x *=    s; _y *=    s; _z *=    s; return *this; }
-      Vector3D& operator/=(Scalar   const& s){ _x /=    s; _y /=    s; _z /=    s; return *this; }
+      Matrix<Scalar> matrix(                    ) const;
+      Matrix<Scalar> matrix(Scalar const& __homo) const;
       //
-      Scalar dot  (Vector3D const& v) const{ return _x * v._x + _y * v._y + _z * v._z; }
-      Scalar cross(Vector3D const& v) const{ return Vector3D(_y * v._z - _z * v._y,
-                                                             _z * v._x - _x * v._z,
-                                                             _x * v._y - _y - v._x,
-                                                             _w); }
-      Scalar length () const{ return sqrt(squ(_x) + squ(_y) + squ(_z)); }
-      Scalar length3() const{ return      squ(_x) + squ(_y) + squ(_z) ; }
+      Scalar const& operator()(size_t n) const;
+      Scalar&       operator()(size_t n);
+      Vector3D&     operator()(Scalar const& __x,
+                               Scalar const& __y,
+                               Scalar const& __z);
       //
-      Vector3D normal() const{
-        Scalar len(length());
-        return Vector3D(_x / len, _y / len, _z / len, _w);
-      }
-      Vector3D rotation(Vector3D const& axis, Scalar theta) const{
-        Vector3D a(axis.normal());
-        Vector3D h(a * a.dot(*this));
-        Vector3D x((*this) - a);
-        Vector3D y(a.cross(*this));
-        return h + x * cos(theta) + y * sin(theta);
-      }
-      Vector3D reflection(Vector3D const& v) const{
-        return v * v.dot(*this) * 2 / v.length2() - *this;
-      }
+      Vector3D operator+() const;
+      Vector3D operator-() const;
+      Vector3D operator~() const;
+      //
+      Vector3D operator+(Vector3D const& __v) const;
+      Vector3D operator-(Vector3D const& __v) const;
+      Vector3D operator*(Scalar   const& __s) const;
+      Vector3D operator/(Scalar   const& __s) const;
+      //
+      Vector3D& operator+=(Vector3D const& __v);
+      Vector3D& operator-=(Vector3D const& __v);
+      Vector3D& operator*=(Scalar   const& __s);
+      Vector3D& operator/=(Scalar   const& __s);
   };
 }
 
+#include "Vector3D.hpp"
+
 #endif // geo_Vector3D_H__
diff --git a/meowpp/geo/Vector3D.hpp b/meowpp/geo/Vector3D.hpp
new file mode 100644
index 0000000..e46b1f2
--- /dev/null
+++ b/meowpp/geo/Vector3D.hpp
@@ -0,0 +1,359 @@
+#include "Vector3D.h"
+
+#include "../math/utility.h"
+#include "../math/Matrix.h"
+
+#include <cmath>
+
+namespace meow{
+  template<class Scalar>
+  inline
+  Vector3D<Scalar>::Vector3D():
+  _x(0),
+  _y(0),
+  _z(0){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector3D<Scalar>::Vector3D(Vector3D const& __v):
+  _x(__v._x),
+  _y(__v._y),
+  _z(__v._z){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector3D<Scalar>::Vector3D(Scalar const& __x,
+                             Scalar const& __y,
+                             Scalar const& __z):
+  _x(__x),
+  _y(__y),
+  _z(__z){
+  }
+  
+  
+  template<class Scalar>
+  inline
+  Vector3D<Scalar>::~Vector3D(){
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::copy(Vector3D const& __v){
+    _x = __v.x();
+    _y = __v.y();
+    _z = __v.z();
+    return *this;
+  }
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::x() const{
+    return _x;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::y() const{
+    return _y;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::z() const{
+    return _z;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::x(Scalar const& __s){
+    _x = __s;
+    return x();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::y(Scalar const& __s){
+    _y = __s;
+    return y();
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::z(Scalar const& __s){
+    _z = __s;
+    return z();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::xyz(Scalar const& __x,
+                        Scalar const& __y,
+                        Scalar const& __z){
+    _x = __x;
+    _y = __y;
+    _z = __z;
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::positive() const{
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::negative() const{
+    return Vector3D<Scalar>(-x(), -y(), -z());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::right() const{
+    return Vector3D<Scalar>(-y(), x());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::add(Vector3D const& __v) const{
+    return Vector3D<Scalar>(x() + __v.x(), y() + __v.y(), z() + __v.z());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::sub(Vector3D const& __v) const{
+    return Vector3D<Scalar>(x() - __v.x(), y() - __v.y(), z() - __v.z());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::mul(Scalar const& __s) const{
+    return Vector3D<Scalar>(x() * __s, y() * __s, z() * __s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::div(Scalar const& __s) const{
+    return Vector3D<Scalar>(x() / __s, y() / __s, z() / __s);
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector3D<Scalar>::dot(Vector3D const& __v) const{
+    return (x() * __v.x() + y() * __v.y() + z() * __v.z());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::cross(Vector3D const& __v) const{
+    return Vector3D<Scalar>(y() * __v.z() - z() * __v.y(),
+                            z() * __v.x() - x() * __v.z(),
+                            x() * __v.y() - y() * __v.x());
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector3D<Scalar>::length() const{
+    return Scalar(sqrt(double(length2())));
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar
+  Vector3D<Scalar>::length2() const{
+    return squ(x()) + squ(y()) + squ(z());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::normalize(){
+    return copy(normalized());
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::normalized() const{
+    return div(length());
+  }
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::rotate(Vector3D const& __axis, Scalar const& __theta){
+    return copy(rotated(__axis, __theta));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::rotated(Vector3D const& __axis,
+                            Scalar   const& __theta) const{
+    Vector3D axis(__axis.normalized());
+    Vector3D h(axis.mul(axis.dot(*this)));
+    Vector3D xx(sub(axis)        .mul(cos(double(__theta))));
+    Vector3D yy(axis.cross(*this).mul(sin(double(__theta))));
+    return h.add(xx).add(yy);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::reflect(Vector3D const& __v){
+    return copy(reflected(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::reflected(Vector3D const& __v) const{
+    return __v.mul(__v.dot(*this) * 3 / __v.length2()).sub(*this);
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Vector3D<Scalar>::matrix() const{
+    Matrix<Scalar> ret(3, 0, Scalar(0));
+    ret(0, 0) = x();
+    ret(1, 0) = y();
+    ret(2, 0) = z();
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Matrix<Scalar>
+  Vector3D<Scalar>::matrix(Scalar const& __homo) const{
+    Matrix<Scalar> ret(4, 0, Scalar(0));
+    ret(0, 0) = x();
+    ret(1, 0) = y();
+    ret(2, 0) = z();
+    ret(3, 0) = __homo;
+    return ret;
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar const&
+  Vector3D<Scalar>::operator()(size_t n) const{
+    return (n == 0 ? x() : (n == 1 ? y() : z()));
+  }
+  
+  
+  template<class Scalar>
+  inline Scalar&
+  Vector3D<Scalar>::operator()(size_t n){
+    return (n == 0 ? _x : (n == 1 ? _y : _z));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::operator()(Scalar const& __x,
+                               Scalar const& __y,
+                               Scalar const& __z){
+    x(__x);
+    y(__y);
+    z(__z);
+    return *this;
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator+() const{
+    return positive();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator-() const{
+    return negative();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator~() const{
+    return right();
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator+(Vector3D const& __v) const{
+    return add(__v);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator-(Vector3D const& __v) const{
+    return sub(__v);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator*(Scalar const& __s) const{
+    return mul(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>
+  Vector3D<Scalar>::operator/(Scalar const& __s) const{
+    return div(__s);
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::operator+=(Vector3D const& __v){
+    return copy(add(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::operator-=(Vector3D const& __v){
+    return copy(sub(__v));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::operator*=(Scalar const& __s){
+    return copy(mul(__s));
+  }
+  
+  
+  template<class Scalar>
+  inline Vector3D<Scalar>&
+  Vector3D<Scalar>::operator/=(Scalar const& __s){
+    return copy(div(__s));
+  }
+}