///////////////////////////////////////////////////////////////////////////////
//
//      ##  ######          
//       ######  ###        Based on:
//  ## ###############      
//   ########## # # #       Shark 3D Engine (www.shark3d.com)
//    ########
//   ######### # # #        Copyright (c) 1996-2001 Spinor GmbH.
//  ##   ##########         
//      ##                  
//        
///////////////////////////////////////////////////////////////////////////////
//
// G A M E.O.N.E - LOW LEVEL LIB V1.0
// Copyright (C) 2001 LEVEL ONE ENTERTAINMENT, 
// Licensed under the terms of LGPL.
//:---------------------------------------------------------------------------
//:Description
//
// LOW LEVEL 3D VEC4 UTIL HEADERFILE
//
//:---------------------------------------------------------------------------

#ifndef _LL3D_VEC4_H
#define _LL3D_VEC4_H


#include "ll3d_math.h"

namespace ll3d {


template<typename T> struct CVec4;
template<typename T> struct CMat4x4;


typedef CVec4<float>	 CVec4f;
typedef const CVec4f	 CVec4f_c;
typedef CVec4f			&CVec4f_r;
typedef const CVec4f	&CVec4f_cr;
typedef CVec4<double>	 CVec4d;
typedef const CVec4d	 CVec4d_c;
typedef CVec4d			&CVec4d_r;
typedef const CVec4d	&CVec4d_cr;

typedef CMat4x4<float>	 CMat4x4f;
typedef const CMat4x4f	 CMat4x4f_c;
typedef CMat4x4f		&CMat4x4f_r;
typedef const CMat4x4f	&CMat4x4f_cr;
typedef CMat4x4<double>  CMat4x4d;
typedef const CMat4x4d	 CMat4x4d_c;
typedef CMat4x4d		&CMat4x4d_r;
typedef const CMat4x4d	&CMat4x4d_cr;

/*
 * 4D Vector
 */

template<typename T>
struct CVec4
{
    T m_x, m_y, m_z, m_w;

    CVec4<T>();

    // Set each component to the same value.
    // Also allows to use 0 as null vector 
    // in the same way as for matrices CMat4x4.
    CVec4<T>(T s);

    // Constructor.
    CVec4<T>(T x, T y, T z, T w);
    CVec4<T>(const CVec3<T> &v, T w);

    // Copy constructor.
    CVec4<T>(const CVec4<T> &v);

    // Assignment operator.
    CVec4<T> &operator=(const CVec4<T> &v);

    CVec4<T> operator*(const CMat4x4<T> &m) const;

    
    bool operator==(const CVec4<T> &v) const;// Equality.

    bool operator!=(const CVec4<T> &v) const;// Inequality.

    
    CVec3<T> GetSubV() const;// Return the 3d sub-vector (x,y,z).
    float GetSubW() const;// Return the w-component.

    void SetSubV(CVec3f_cr SubV);// Set the sub-vector (x,y,z).
    void SetSubW(float SubW);// Set the w-component.

    CVec3<T> GetProj() const;// Return the projected vector (x,y,z) / w.
};

/*
 * 4x4 Matrix
 */

template<typename T>
struct CMat4x4
{
    T m_xx, m_xy, m_xz, m_xw,
      m_yx, m_yy, m_yz, m_yw,
      m_zx, m_zy, m_zz, m_zw,
      m_wx, m_wy, m_wz, m_ww;

    CMat4x4<T>();

    // Initializating with multiple of unity matrix.
    // Also allows to use 0 as null matrix and 1 as unity matrix.
    CMat4x4<T>(T s);

    // Constructor.
    CMat4x4<T>(
        T xx, T xy, T xz, T xw,
        T yx, T yy, T yz, T yw,
        T zx, T zy, T zz, T zw,
        T wx, T wy, T wz, T ww);

    // Construct by the matrix CMat4x4::GetSubVV,
    // the vectors CMat4x4::GetSubVW 
    // and CMat4x4::GetSubWV, and CMat4x4::GetSubWW.
    CMat4x4<T>(const CMat3x3<T> &vv, const CVec3<T> &vw,
                    const CVec3<T> &wv, T ww);

    
    CMat4x4<T>(const CMat4x4<T> &m);// Copy constructor.
    CMat4x4<T> &operator=(const CMat4x4<T> &m);// Assignment operator.
    CMat4x4<T> operator+(const CMat4x4<T> &m) const;// Matrix addition.
    CMat4x4<T> operator-(const CMat4x4<T> &m) const;// Matrix substraction.
    CMat4x4<T> operator*(const CMat4x4<T> &m) const;// Matrix multiplication.
    CVec4<T> operator*(const CVec4<T> &v) const;
    bool operator==(const CMat4x4<T> &m) const;// Equality.
    bool operator!=(const CMat4x4<T> &m) const;// Inequality.

    CMat4x4<T> GetTransposed() const;// Returns the transposed matrix.

    // Calculate a 3x3 determinant.
    static T Get3x3Det(T xx, T xy, T xz, T yx, T yy, T yz, T zx, T zy, T zz);

    
    CMat4x4<T> GetInv() const;// Return the inverse matrix.

    
    T GetDet() const;// Returns the determinant.

    // Return sub-matrix ((xx, xy, xz), (yx, yy, yz), (zx, zy, zz)).
    CMat3x3<T> GetSubVV() const;

    
    CVec3<T> GetSubVW() const;// Return sub-vector (xw, yw, zw).
    CVec3<T> GetSubWV() const;// Return sub-vector (wx, wy, wz).
    T GetSubWW() const;// Return component ww.

    // Set the sub-matrix ((xx, xy, xz), (yx, yy, yz), (zx, zy, zz)).
    void SetSubVV(const CMat3x3<T> &SubVV);

    
    void SetSubVW(const CVec3<T> &SubVW);// Set the sub-vector (xw, yw, zw).
    void SetSubWV(const CVec3<T> &SubWV);// Set the sub-vector (wx, wy, wz).

    void SetSubWW(T SubWW);// Set the component ww.
    CVec3<T> ProdOneV(const CVec3<T> &v) const;// m.ProdOneV(v) = (m * CVec4(v, 1)).GetSubV();

    T ProdOneW(const CVec3<T> &v) const;// m.ProdOneW(v) = (m * CVec4(v, 1)).m_w;

    
    CVec4<T> ProdOne(const CVec3<T> &v) const;// m.ProdOne(v) = m * CVec4(v, 1);

    // m.ProjOne(v) 
    // = m.@sCMat4x4::ProdOne(v).@sCVec4::GetProj();
    CVec3<T> ProjOne(const CVec3<T> &v) const;

    // Derivation of CMat4x4::ProjOne at v.
    CMat3x3<T> ProjOneDeriv(const CVec3<T> &v) const;

    // m.ProdZeroV(v) = (m * CVec4(v, 0)).GetSubV();
    CVec3<T> ProdZeroV(const CVec3<T> &v) const;

    // m.ProdZeroW(v) = (m * CVec4(v, 0)).m_w;
    T ProdZeroW(const CVec3<T> &v) const;

    // m.ProdZero(v) = m * CVec4(v, 0))
    CVec4<T> ProdZero(const CVec3<T> &v) const;
};

/*
 * Implementation
 */

template<typename T>
inline CVec4<T>::CVec4()
{
}

template<typename T>
inline CVec4<T>::CVec4(T s)
{
    m_x = m_y = m_z = m_w = s;
}

template<typename T>
inline CVec4<T>::CVec4(T x, T y, T z, T w)
{
    m_x = x;
    m_y = y;
    m_z = z;
    m_w = w;
}

template<typename T>
CVec4<T>::CVec4(const CVec3<T> &v, T w)
{
    m_x = v.m_x;
    m_y = v.m_y;
    m_z = v.m_z;
    m_w = w;
}

template<typename T>
inline CVec4<T>::CVec4(const CVec4<T> &v)
{
    m_x = v.m_x;
    m_y = v.m_y;
    m_z = v.m_z;
    m_w = v.m_w;
}

template<typename T>
inline CVec4<T> &CVec4<T>::operator=(const CVec4<T> &v)
{
    m_x = v.m_x;
    m_y = v.m_y;
    m_z = v.m_z;
    m_w = v.m_w;
    return *this;
}

template<typename T>
CVec4<T> CVec4<T>::operator*(const CMat4x4<T> &m) const
{
    return CVec4<T>(
        m_x * m.m_xx + m_y * m.m_yx + m_z * m.m_zx + m_w * m.m_wx,
        m_x * m.m_xy + m_y * m.m_yy + m_z * m.m_zy + m_w * m.m_wy,
        m_x * m.m_xz + m_y * m.m_yz + m_z * m.m_zz + m_w * m.m_wz,
        m_x * m.m_xw + m_y * m.m_yw + m_z * m.m_zw + m_w * m.m_ww);
}

template<typename T>
inline bool CVec4<T>::operator==(const CVec4<T> &v) const
{
    return m_x == v.m_x && m_y == v.m_y && m_z == v.m_z && m_w == v.m_w;
}

template<typename T>
inline bool CVec4<T>::operator!=(const CVec4<T> &v) const
{
    return m_x != v.m_x || m_y != v.m_y || m_z != v.m_z || m_w != v.m_w;
}

template<typename T>
CVec3<T> CVec4<T>::GetSubV() const
{
    return CVec3<T>(m_x, m_y, m_z);
}

template<typename T>
float CVec4<T>::GetSubW() const
{
    return m_w;
}

template<typename T>
void CVec4<T>::SetSubV(CVec3f_cr SubV)
{
    m_x = SubV.m_x;
    m_y = SubV.m_y;
    m_z = SubV.m_z;
}

template<typename T>
void CVec4<T>::SetSubW(float SubW)
{
    m_w = SubW;
}

template<typename T>
CVec3<T> CVec4<T>::GetProj() const
{
    return CVec3<T>(m_x, m_y, m_z) / m_w;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
inline CMat4x4<T>::CMat4x4()
{
}

template<typename T>
inline CMat4x4<T>::CMat4x4(T s)
{
    m_xx = s;   m_xy = 0;   m_xz = 0;   m_xw = 0;
    m_yx = 0;   m_yy = s;   m_yz = 0;   m_yw = 0;
    m_zx = 0,   m_zy = 0;   m_zz = s;   m_zw = 0;
    m_wx = 0,   m_wy = 0;   m_wz = 0;   m_ww = s;
}

template<typename T>
inline CMat4x4<T>::CMat4x4(T xx, T xy, T xz, T xw,
                                      T yx, T yy, T yz, T yw,
                                      T zx, T zy, T zz, T zw,
                                      T wx, T wy, T wz, T ww)
{
    m_xx = xx;  m_xy = xy;  m_xz = xz;  m_xw = xw;
    m_yx = yx;  m_yy = yy;  m_yz = yz;  m_yw = yw;
    m_zx = zx;  m_zy = zy;  m_zz = zz;  m_zw = zw;
    m_wx = wx;  m_wy = wy;  m_wz = wz;  m_ww = ww;
}

template<typename T>
CMat4x4<T>::CMat4x4(
                    const CMat3x3<T> &vv, const CVec3<T> &vw,
                    const CVec3<T> &wv, T ww)
{
    m_xx = vv.m_xx; m_xy = vv.m_xy; m_xz = vv.m_xz; m_xw = vw.m_x;
    m_yx = vv.m_yx; m_yy = vv.m_yy; m_yz = vv.m_yz; m_yw = vw.m_y;
    m_zx = vv.m_zx; m_zy = vv.m_zy; m_zz = vv.m_zz; m_zw = vw.m_z;
    m_wx = wv.m_x;  m_wy = wv.m_y;  m_wz = wv.m_z;  m_ww = ww;
}

template<typename T>
inline CMat4x4<T>::CMat4x4(const CMat4x4<T> &m)
{
    m_xx = m.m_xx;  m_xy = m.m_xy;  m_xz = m.m_xz;  m_xw = m.m_xw;
    m_yx = m.m_yx;  m_yy = m.m_yy;  m_yz = m.m_yz;  m_yw = m.m_yw;
    m_zx = m.m_zx;  m_zy = m.m_zy;  m_zz = m.m_zz;  m_zw = m.m_zw;
    m_wx = m.m_wx;  m_wy = m.m_wy;  m_wz = m.m_wz;  m_ww = m.m_ww;
}

template<typename T>
inline CMat4x4<T> &CMat4x4<T>::operator=(const CMat4x4<T> &m)
{
    m_xx = m.m_xx;  m_xy = m.m_xy;  m_xz = m.m_xz;  m_xw = m.m_xw;
    m_yx = m.m_yx;  m_yy = m.m_yy;  m_yz = m.m_yz;  m_yw = m.m_yw;
    m_zx = m.m_zx;  m_zy = m.m_zy;  m_zz = m.m_zz;  m_zw = m.m_zw;
    m_wx = m.m_wx;  m_wy = m.m_wy;  m_wz = m.m_wz;  m_ww = m.m_ww;
    return *this;
}

///////////////////////////////////////////////////////////////////////////////

template<typename T>
CMat4x4<T> CMat4x4<T>::operator+(const CMat4x4<T> &m) const
{
    return CMat4x4<T>(
        m_xx + m.m_xx, m_xy + m.m_xy, m_xz + m.m_xz, m_xw + m.m_xw,
        m_yx + m.m_yx, m_yy + m.m_yy, m_yz + m.m_yz, m_xw + m.m_yw,
        m_zx + m.m_zx, m_zy + m.m_zy, m_zz + m.m_zz, m_xw + m.m_zw,
        m_wx + m.m_wx, m_wy + m.m_wy, m_wz + m.m_wz, m_xw + m.m_ww);
}

template<typename T>
CMat4x4<T> CMat4x4<T>::operator-(const CMat4x4<T> &m) const
{
    return CMat4x4<T>(
        m_xx - m.m_xx, m_xy - m.m_xy, m_xz - m.m_xz, m_xw - m.m_xw,
        m_yx - m.m_yx, m_yy - m.m_yy, m_yz - m.m_yz, m_xw - m.m_yw,
        m_zx - m.m_zx, m_zy - m.m_zy, m_zz - m.m_zz, m_xw - m.m_zw,
        m_wx - m.m_wx, m_wy - m.m_wy, m_wz - m.m_wz, m_xw - m.m_ww);
}

template<typename T>
CMat4x4<T> CMat4x4<T>::operator*(const CMat4x4<T> &m) const
{
    return CMat4x4<T>(
        m_xx * m.m_xx + m_xy * m.m_yx + m_xz * m.m_zx + m_xw * m.m_wx,
        m_xx * m.m_xy + m_xy * m.m_yy + m_xz * m.m_zy + m_xw * m.m_wy,
        m_xx * m.m_xz + m_xy * m.m_yz + m_xz * m.m_zz + m_xw * m.m_wz,
        m_xx * m.m_xw + m_xy * m.m_yw + m_xz * m.m_zw + m_xw * m.m_ww,
        m_yx * m.m_xx + m_yy * m.m_yx + m_yz * m.m_zx + m_yw * m.m_wx,
        m_yx * m.m_xy + m_yy * m.m_yy + m_yz * m.m_zy + m_yw * m.m_wy,
        m_yx * m.m_xz + m_yy * m.m_yz + m_yz * m.m_zz + m_yw * m.m_wz,
        m_yx * m.m_xw + m_yy * m.m_yw + m_yz * m.m_zw + m_yw * m.m_ww,
        m_zx * m.m_xx + m_zy * m.m_yx + m_zz * m.m_zx + m_zw * m.m_wx,
        m_zx * m.m_xy + m_zy * m.m_yy + m_zz * m.m_zy + m_zw * m.m_wy,
        m_zx * m.m_xz + m_zy * m.m_yz + m_zz * m.m_zz + m_zw * m.m_wz,
        m_zx * m.m_xw + m_zy * m.m_yw + m_zz * m.m_zw + m_zw * m.m_ww,
        m_wx * m.m_xx + m_wy * m.m_yx + m_wz * m.m_zx + m_ww * m.m_wx,
        m_wx * m.m_xy + m_wy * m.m_yy + m_wz * m.m_zy + m_ww * m.m_wy,
        m_wx * m.m_xz + m_wy * m.m_yz + m_wz * m.m_zz + m_ww * m.m_wz,
        m_wx * m.m_xw + m_wy * m.m_yw + m_wz * m.m_zw + m_ww * m.m_ww);
}

template<typename T>
CVec4<T> CMat4x4<T>::operator*(const CVec4<T> &m) const
{
    return CVec4<T>(
        m_xx * m.m_x + m_xy * m.m_y + m_xz * m.m_z + m_xw * m.m_w,
        m_yx * m.m_x + m_yy * m.m_y + m_yz * m.m_z + m_yw * m.m_w,
        m_zx * m.m_x + m_zy * m.m_y + m_zz * m.m_z + m_zw * m.m_w,
        m_wx * m.m_x + m_wy * m.m_y + m_wz * m.m_z + m_ww * m.m_w);
}

template<typename T>
bool CMat4x4<T>::operator==(const CMat4x4<T> &m) const
{
    return 
        m_xx == m.m_xx && m_xy == m.m_xy && m_xz == m.m_xz && m_xw == m.m_xw
     && m_yx == m.m_yx && m_yy == m.m_yy && m_yz == m.m_yz && m_yw == m.m_yw
     && m_zx == m.m_zx && m_zy == m.m_zy && m_zz == m.m_zz && m_zw == m.m_zw
     && m_wx == m.m_wx && m_wy == m.m_wy && m_wz == m.m_wz && m_ww == m.m_ww;
}

template<typename T>
bool CMat4x4<T>::operator!=(const CMat4x4<T> &m) const
{
    return 
        m_xx != m.m_xx || m_xy != m.m_xy || m_xz != m.m_xz || m_xw != m.m_xw
     || m_yx != m.m_yx || m_yy != m.m_yy || m_yz != m.m_yz || m_yw != m.m_yw
     || m_zx != m.m_zx || m_zy != m.m_zy || m_zz != m.m_zz || m_zw != m.m_zw
     || m_wx != m.m_wx || m_wy != m.m_wy || m_wz != m.m_wz || m_ww != m.m_ww;
}

template<typename T>
T CMat4x4<T>::Get3x3Det(
                T xx, T xy, T xz, T yx, T yy, T yz, T zx, T zy, T zz)
{
    return xx * yy * zz - xx * yz * zy 
         + xy * yz * zx - xy * yx * zz 
         + xz * yx * zy - xz * yy * zx;
}

template<typename T>
inline CMat4x4<T> CMat4x4<T>::GetTransposed() const
{
    return CMat4x4<T>(
        m_xx, m_yx, m_zx, m_wx,
        m_xy, m_yy, m_zy, m_wy,
        m_xz, m_yz, m_zz, m_wz,
        m_xw, m_yw, m_zw, m_ww);
}

template<typename T>
CMat4x4<T> CMat4x4<T>::GetInv() const
{
    T sxx = + Get3x3Det(m_yy, m_yz, m_yw, m_zy, m_zz, m_zw, m_wy, m_wz, m_ww);
    T sxy = - Get3x3Det(m_yz, m_yw, m_yx, m_zz, m_zw, m_zx, m_wz, m_ww, m_wx);
    T sxz = + Get3x3Det(m_yw, m_yx, m_yy, m_zw, m_zx, m_zy, m_ww, m_wx, m_wy);
    T sxw = - Get3x3Det(m_yx, m_yy, m_yz, m_zx, m_zy, m_zz, m_wx, m_wy, m_wz);
    T syx = - Get3x3Det(m_zy, m_zz, m_zw, m_wy, m_wz, m_ww, m_xy, m_xz, m_xw);
    T syy = + Get3x3Det(m_zz, m_zw, m_zx, m_wz, m_ww, m_wx, m_xz, m_xw, m_xx);
    T syz = - Get3x3Det(m_zw, m_zx, m_zy, m_ww, m_wx, m_wy, m_xw, m_xx, m_xy);
    T syw = + Get3x3Det(m_zx, m_zy, m_zz, m_wx, m_wy, m_wz, m_xx, m_xy, m_xz);
    T szx = + Get3x3Det(m_wy, m_wz, m_ww, m_xy, m_xz, m_xw, m_yy, m_yz, m_yw);
    T szy = - Get3x3Det(m_wz, m_ww, m_wx, m_xz, m_xw, m_xx, m_yz, m_yw, m_yx);
    T szz = + Get3x3Det(m_ww, m_wx, m_wy, m_xw, m_xx, m_xy, m_yw, m_yx, m_yy);
    T szw = - Get3x3Det(m_wx, m_wy, m_wz, m_xx, m_xy, m_xz, m_yx, m_yy, m_yz);
    T swx = - Get3x3Det(m_xy, m_xz, m_xw, m_yy, m_yz, m_yw, m_zy, m_zz, m_zw);
    T swy = + Get3x3Det(m_xz, m_xw, m_xx, m_yz, m_yw, m_yx, m_zz, m_zw, m_zx);
    T swz = - Get3x3Det(m_xw, m_xx, m_xy, m_yw, m_yx, m_yy, m_zw, m_zx, m_zy);
    T sww = + Get3x3Det(m_xx, m_xy, m_xz, m_yx, m_yy, m_yz, m_zx, m_zy, m_zz);
    T d = m_xx * sxx + m_xy * sxy + m_xz * sxz + m_xw * sxw;
    T f = 1 / d;
    return CMat4x4<T>(
        f * sxx, f * syx, f * szx, f * swx,
        f * sxy, f * syy, f * szy, f * swy,
        f * sxz, f * syz, f * szz, f * swz,
        f * sxw, f * syw, f * szw, f * sww);
}

template<typename T>
T CMat4x4<T>::GetDet() const
{
    T sxx = Get3x3Det(m_yy, m_yz, m_yw, m_zy, m_zz, m_zw, m_wy, m_wz, m_ww);
    T sxy = Get3x3Det(m_yz, m_yw, m_yx, m_zz, m_zw, m_zx, m_wz, m_ww, m_wx);
    T sxz = Get3x3Det(m_yw, m_yx, m_yy, m_zw, m_zx, m_zy, m_ww, m_wx, m_wy);
    T sxw = Get3x3Det(m_yx, m_yy, m_yz, m_zx, m_zy, m_zz, m_wx, m_wy, m_wz);
    return m_xx * sxx + m_xy * sxy + m_xz * sxz + m_xw * sxw;
}

template<typename T>
CMat3x3<T> CMat4x4<T>::GetSubVV() const
{
    return CMat3x3<T>(
                m_xx, m_xy, m_xz, m_yx, m_yy, m_yz, m_zx, m_zy, m_zz);
}

template<typename T>
CVec3<T> CMat4x4<T>::GetSubVW() const
{
    return CVec3<T>(m_xw, m_yw, m_zw);
}

template<typename T>
CVec3<T> CMat4x4<T>::GetSubWV() const
{
    return CVec3<T>(m_wx, m_wy, m_wz);
}

template<typename T>
T CMat4x4<T>::GetSubWW() const
{
    return m_ww;
}

template<typename T>
void CMat4x4<T>::SetSubVV(const CMat3x3<T> &SubVV)
{
    m_xx = SubVV.m_xx; m_xy = SubVV.m_xy; m_xz = SubVV.m_xz;
    m_yx = SubVV.m_yx; m_yy = SubVV.m_yy; m_yz = SubVV.m_yz;
    m_zx = SubVV.m_zx; m_zy = SubVV.m_zy; m_zz = SubVV.m_zz;
}

template<typename T>
void CMat4x4<T>::SetSubVW(const CVec3<T> &SubVW)
{
    m_xw = SubVW.m_x;
    m_yw = SubVW.m_y;
    m_zw = SubVW.m_z;
}

template<typename T>
void CMat4x4<T>::SetSubWV(const CVec3<T> &SubWV)
{
    m_wx = SubWV.m_x;
    m_wy = SubWV.m_y;
    m_wz = SubWV.m_z;
}

template<typename T>
void CMat4x4<T>::SetSubWW(T SubWW)
{
    m_ww = SubWW;
}

template<typename T>
CVec3<T> CMat4x4<T>::ProdOneV(const CVec3<T> &v) const
{
    return CVec3<T>(
        m_xx * v.m_x + m_xy * v.m_y + m_xz * v.m_z + m_xw,
        m_yx * v.m_x + m_yy * v.m_y + m_yz * v.m_z + m_yw,
        m_zx * v.m_x + m_zy * v.m_y + m_zz * v.m_z + m_zw);
}

template<typename T>
T CMat4x4<T>::ProdOneW(const CVec3<T> &v) const
{
    return m_wx * v.m_x + m_wy * v.m_y + m_wz * v.m_z + m_ww;
}

template<typename T>
CVec4<T> CMat4x4<T>::ProdOne(const CVec3<T> &v) const
{
    return CVec4<T>(
        m_xx * v.m_x + m_xy * v.m_y + m_xz * v.m_z + m_xw,
        m_yx * v.m_x + m_yy * v.m_y + m_yz * v.m_z + m_yw,
        m_zx * v.m_x + m_zy * v.m_y + m_zz * v.m_z + m_zw,
        m_wx * v.m_x + m_wy * v.m_y + m_wz * v.m_z + m_ww);
}

template<typename T>
CVec3<T> CMat4x4<T>::ProjOne(const CVec3<T> &v) const
{
    return ProdOne(v).GetProj();
}

template<typename T>
CMat3x3<T> CMat4x4<T>::ProjOneDeriv(const CVec3<T> &v) const
{
    CVec3<T> LinV = ProdOneV(v);
    T LinW = ProdOneW(v);

    CMat3x3<T> Num = GetSubVV() * LinW - CMat3x3<T>::ByTensorProd(LinV, GetSubWV());
    return Num / (LinW * LinW);
}

template<typename T>
CVec3<T> CMat4x4<T>::ProdZeroV(const CVec3<T> &v) const
{
    return CVec3<T>(
        m_xx * v.m_x + m_xy * v.m_y + m_xz * v.m_z,
        m_yx * v.m_x + m_yy * v.m_y + m_yz * v.m_z,
        m_zx * v.m_x + m_zy * v.m_y + m_zz * v.m_z);
}

template<typename T>
T CMat4x4<T>::ProdZeroW(const CVec3<T> &v) const
{
    return m_wx * v.m_x + m_wy * v.m_y + m_wz * v.m_z;
}

template<typename T>
CVec4<T> CMat4x4<T>::ProdZero(const CVec3<T> &v) const
{
    return CVec4<T>(
        m_xx * v.m_x + m_xy * v.m_y + m_xz * v.m_z,
        m_yx * v.m_x + m_yy * v.m_y + m_yz * v.m_z,
        m_zx * v.m_x + m_zy * v.m_y + m_zz * v.m_z,
        m_wx * v.m_x + m_wy * v.m_y + m_wz * v.m_z);
}


}//namespace ll3d
#endif //_LL3D_VEC3_H