(Due to a grievance filed by a student in one of my classes, the following disclaimer has been patterned from a Microsoft web page and is applied to code found herein on Dr. Duchowski's web pages.)

WARNING: ANY USE BY YOU OF THE CODE PROVIDED IN THESE WEB PAGES IS AT YOUR OWN RISK. Dr. Duchowski provides this code "as is" without warranty of any kind, either express or implied, including but not limited to the implied warranties of merchantability and/or fitness for a particular purpose.

Dr. Duchowski will not be held responsible for any inaccuracies, bugs, documented or undocumented, in the sample code below. As originally stated in the first line of this web page, these are merely suggestions. SHOULD YOU WISH TO USE THESE EXAMPLES FOR COURSEWORK, YOU ARE DOING SO AT YOUR OWN RISK. Dr. Duchowski will not grant extra credit, time, debugging time, or anything else should you find a bug in the code (except maybe a "Thanks for pointing that out!").

Original page starts from line below.

Suggestions for implementation of a math library (e.g., mathlib)

Objectives:
To implement a C++ math library suitable for use in a graphics program.

Description:
The main classes are matrix, quaternion, both rely on and use the STL's std::vector<> class. To facilitate various "inter-class" operations (such vector-matrix multiplication), and to promote commonality of operations, the classes are each others' friends.

Your math library should run something like this driver.cpp program, generating mathematically correct results.

When testing your math class, you should probably start with the vector class, redefining some of the operators of the class (e.g., operator<<, operator- (both unary and binary), operator*, operator+, and define member functions for the dot and cross products, and normalization.

The matrix class:
Define the matrix class so that it can be used to generally represent matrices of any size (although we mainly need 4x4 matrices). Since OpenGL can accept pointers to arrays, you can create the class with a private data member consisting of a simple array. Alternatively, you can use an array of std::vector<> like this: 
template 
class matrix {
	...
        private:
        std::vector<std::vector<T> >    arr;
};
If you use this, you can then use the operator[] to get at matrix rows and individual elements via this piece of code I got from Mark A. Weiss' book (Data Structures and Algorithm Analysis in C++, 3rd ed., Pearson Education (Addison Wesley), Boston, MA, 2006): 
        // operators (dereferencing)
        const std::vector<T> & operator[](int r) const  { return arr[r]; }
        std::vector<T> & operator[](int r)              { return arr[r]; }
};

Note that the above matrix representation may not be very useful for OpenGL single-array parameter passing, however I like it for its reliance on the STL std::vector<>. It seemed elegant, and it allows using all of the functions that are already provided for that class (e.g., size() etc.).

The following are recommended member functions for the matrix (most are self-explanatory): 

        public:
        matrix(int r=4, int c=4) : arr(r)
          { for(int i=0;i<r;++i) arr[i].resize(c); }
        matrix(const matrix& rhs) : arr(rhs.rows())
          { for(int i=0;i<rows();++i) arr[i] = rhs.arr[i]; }

        // destructors – default should be ok

        // operators (dereferencing)
        const std::vector<T> & operator[](int r) const  { return arr[r]; }
        std::vector<T> & operator[](int r)              { return arr[r]; }

        // operators (copy assignment)
        matrix  operator=(const matrix& rhs);

        // operators (unary)
        matrix  operator-(void);
        matrix  operator!(void);

        // operators (const matrix&, const matrix&)
        int     operator==(const matrix& rhs);
        int     operator<(const matrix& rhs);
        int     operator>(const matrix& rhs);
        matrix  operator*=(const matrix& rhs);
        matrix  operator*(const matrix& rhs);
        matrix  operator/=(matrix& rhs);
        matrix  operator/(matrix& rhs);
        matrix  operator+=(const matrix& rhs);
        matrix  operator+(const matrix& rhs);
        matrix  operator-=(const matrix& rhs);
        matrix  operator-(const matrix& rhs);

        // operators (const matrix&, scalar)
        matrix  operator*=(const T& rhs);
        matrix  operator*(const T& rhs);
        matrix  operator/=(const T& rhs);
        matrix  operator/(const T& rhs);
        matrix  operator+=(const T& rhs);
        matrix  operator+(const T& rhs);
        matrix  operator-=(const T& rhs);
        matrix  operator-(const T& rhs);

        // operators (const scalar&, const matrix&)
        // note: g++ outputs error 'declaration of `operator*' as non-function'
        // the problem is that the declaration of member operator* hides all
        // overloaded declarations from global scope—to access the global
        // function, you need to explicitly qualify it:
        //   friend matrix ::operator* <> (const T&, const matrix&);
        // but this will not compile since :: gets associated with matrix<T>;
        // instead, put parenthesis around the function name
        //   friend matrix (::operator* <>) (const T&, const matrix&);
#ifdef ANM_OSX
        // since we need to use g++3.3 on the Macs...
        friend matrix operator* <>(const T& lhs, const matrix& rhs);
        friend matrix operator/ <>(const T& lhs, const matrix& rhs);
        friend matrix operator+ <>(const T& lhs, const matrix& rhs);
        friend matrix operator- <>(const T& lhs, const matrix& rhs);

        friend std::vector<T> operator* <>(const std::vector<T>&, const matrix&);
#else
        friend matrix (::operator* <>)(const T& lhs, const matrix& rhs);
        friend matrix (::operator/ <>)(const T& lhs, const matrix& rhs);
        friend matrix (::operator+ <>)(const T& lhs, const matrix& rhs);
        friend matrix (::operator- <>)(const T& lhs, const matrix& rhs);

        friend std::vector<T> (::operator* <>)(const std::vector<T>&, const matrix&);
#endif

        // friends – note the extra <> telling the compiler to instantiate
        // a templated version of operator<< – <T> is also legal, i.e.,
        // friend std::ostream& operator<< <T>(std::ostream&, const matrix&);
//      friend std::istream&    operator>>(std::istream&, matrix&);
//      friend std::istream&    operator>>(std::istream&, matrix*)
//        { return(s >> (*rhs)); }
        friend std::ostream&    operator<< <>(std::ostream&, const matrix&);
        friend std::ostream&    operator<<(std::ostream& s, matrix* rhs)
          { return(s << (*rhs)); }

        // member functions
        int     rows(void) const        { return arr.size(); }
        int     cols(void) const        { return rows() ? arr[0].size() : 0; }
        void    zero(void);
        void    one(void);
        void    identity(void);
        matrix  clamp(T lo,T hi);
        matrix  transpose(void);
        T       trace(void) const;
Note that the above class interface is for a templated matrix class. YOU DO NOT NEED TO USE TEMPLATES. You can safely define your matrix to use std::vector<double>.

The vector class:
Using the STL's std::vector<> is nice since it provides robust functionality. All you have to do is customize the class to your liking. Here's an example of how I've used it so far: 
#ifndef VECTOR_H
#define VECTOR_H
#include <iostream>
#include <vector>

template <typename T>
std::ostream&   operator<<(std::ostream& s, const std::vector<T>&);

template <typename T>
std::vector<T>  operator-(const std::vector<T>&);

template <typename T> 
std::vector<T>  operator*(const T&, const std::vector<T>&);

template <typename T>
std::vector<T>  operator+(const std::vector<T>&, const std::vector<T>&);

template <typename T>
std::vector<T>  operator-(const std::vector<T>&, const std::vector<T>&);

template <typename T> 
T               dot(const std::vector<T>&, const std::vector<T>&); 

template <typename T>
std::vector<T>  cross(const std::vector<T>&, const std::vector<T>&);

template <typename T>
std::vector<T>  norm(const std::vector<T>&);
#endif

The quaternion class:
A quaternion is defined by two components: a scalar value (s) and a vector (v), e.g., q = (s, [v1 v2 v3]). A simple way to implement the quaternion class is to once again to just use the STL's std::vector<>: 
class quaternion {
        private:
        std::vector<T>  q;                      // [x,y,z,s]
};
You just have to remember which element s is, is it the 0th or the 3rd. The following are the recommended operations: 
        public:
        // constructors (overloaded)
        quaternion(T sin=0.0) : \
                q(4,0.0)                        { q[3] = sin; }
        quaternion(const quaternion& rhs) : \
                q(rhs.q)                        { };
        quaternion(quaternion *rhs) : \
                q(rhs->q)                       { };
        quaternion(std::vector<T> qin) : \
                q(4,0.0)                        { q[0] = qin[0];
                                                  q[1] = qin[1];
                                                  q[2] = qin[2];
                                                  q[3] = 1.0;
                                                }
        quaternion(matrix<T>& m) : \
                q(4,0.0)                        { q[3] = 1.0; mattoquat(m); }
        quaternion(T sin, std::vector<T> qin) : \
                q(4,0.0)                        { q[0] = qin[0];
                                                  q[1] = qin[1];
                                                  q[2] = qin[2];
                                                  q[3] = sin;
                                                }
        quaternion(T x1,T y1,T x2,T y2) : \
                q(4,0.0)                        { q[3] = 1.0;
                                                  trackball(x1,y1,x2,y2);
                                                }
        // destructors
        //~quaternion()                         { };

        const T& operator[](int r) const        { return q[r]; }
        T& operator[](int r)                    { return q[r]; }

        // operators: unary
        quaternion operator-(void);
        quaternion operator!(void);                     // quaternion inverse

        // operators: (quaternion, const matrix&)
        quaternion operator*(const matrix<T>& rhs);

        // operators: (quaternion, const quaternion&)
        int operator==(const quaternion& rhs);
        int operator<(const quaternion& rhs);
        int operator>(const quaternion& rhs);
        quaternion operator=(const quaternion& rhs);
        quaternion operator*=(const quaternion& rhs);
        quaternion operator*(const quaternion& rhs);
        quaternion operator/=(quaternion& rhs);
        quaternion operator/(quaternion& rhs);
        quaternion operator+=(const quaternion& rhs);
        quaternion operator+(const quaternion& rhs);
        quaternion operator-=(const quaternion& rhs);
        quaternion operator-(const quaternion& rhs);

        // operators: (quaternion, quaternion*)
        int operator==(quaternion* rhs)         { return((*this) == (*rhs)); }
        int operator<(quaternion* rhs)          { return((*this) == (*rhs)); }
        int operator>(quaternion* rhs)          { return((*this) == (*rhs)); }
        quaternion operator=(quaternion* rhs)   { return((*this)  = (*rhs)); }
        quaternion operator*=(quaternion* rhs)  { return((*this) *= (*rhs)); }
        quaternion operator*(quaternion* rhs)   { return((*this) *  (*rhs)); }
        quaternion operator/=(quaternion* rhs)  { return((*this) /= (*rhs)); }
        quaternion operator/(quaternion* rhs)   { return((*this) /  (*rhs)); }
        quaternion operator+=(quaternion* rhs)  { return((*this) += (*rhs)); }
        quaternion operator+(quaternion* rhs)   { return((*this) +  (*rhs)); }
        quaternion operator-=(quaternion* rhs)  { return((*this) -= (*rhs)); }
        quaternion operator-(quaternion* rhs)   { return((*this) -  (*rhs)); }

        // operators (const quaternion&, scalar)
        quaternion operator*=(const T& rhs);
        quaternion operator*(const T& rhs);
        quaternion operator/=(const T& rhs);
        quaternion operator/(const T& rhs);

        // friend operators: (const scalar*, const quaternion&)
#ifdef ANM_OSX
        friend quaternion operator* <>(T lhs,quaternion& rhs);
        friend quaternion operator/ <>(T lhs,quaternion& rhs);
#else
        friend quaternion (::operator* <>)(T lhs,quaternion& rhs);
        friend quaternion (::operator/ <>)(T lhs,quaternion& rhs);
#endif

        // friends
        friend std::ostream& operator<< <>(std::ostream&, const quaternion&);
//      friend std::ostream& operator<<(std::ostream& s, quaternion* rhs)
//        { return(s << (*rhs)); }

        // member functions
        const T&        s()                     { return(q[3]); }
        std::vector<T>  v()                     { std::vector<T> r(3,0.0);
                                                  for(int i=0;i<3;++i)
                                                    r[i] = q[i];
                                                  return(r);
                                                }
        T               dot(quaternion& rhs);
        T               dot(quaternion* rhs)    { return(dot((*rhs))); }
        T               sqlength(void)          { return(dot(this)); }
        T               length(void)            { return(sqrt(dot(this))); }
        quaternion      norm(void)              { return(*this / length()); }
        matrix<T>       L(void);
        matrix<T>       R(void);
        matrix<T>       quattomat(void);        // return (4x4) rot matrix
        void            mattoquat(matrix<T>& m);
        quaternion      slerp2(quaternion rhs,T t);// spherical linear interp.
        quaternion      slerp(quaternion rhs,T t); // spherical linear interp.

        // protected: available to this (base) class and subs. derived classes
        protected:
        void            trackball(T x1,T y1,T x2,T y2);
        T               project_to_sphere(T r,T x,T y);
For a quaternion defined explicitly, e.g., q = (s,v), note a few of the more important operations:
  1. quaternion inverse is defined as
    q-1 = 1 / |q|2 * (s, -v)
    where
    |q|2 = s2 + (v . v),
    i.e., the quaternion's length squared (or its dot procduct).
  2. quaternion multiplication, e.g., q1 * q2 is defined as
    q1 * q2 = (s1 * s2 - v1 . v2, s1 * v2 + s2 * v1 + v1 x v2)
    where the first part before the comma is a scalar and the remaining part is a vector (note the cross product).
  3. the member quattomat is particularly important since it returns a 4x4 matrix that can be fed more or less directly to OpenGL. The matrix returned is obtained as follows (in row-major order): 
      with q = (s, [x y z])
  return (4x4) rotation matrix, in row-major order:
  +-             -+ +-             -+
  | s   z  -y   x | | s   z  -y  -x |
  |-z   s   x   y | |-z   s   x  -y | =
  | y  -x   s   z | | y  -x   s  -z |
  |-x  -y  -z   s | | x   y   z   s |
  +-             -+ +-             -+
  
  +-                                                                  -+
  | (s2+x2-y2-z2)     (2xy+2sz)         (2xz-2sy)             0         |
  |     (2xy-2sz)     (s2-x2+y2-z2)     (2yz+2sx)             0         |
  |     (2xz+2sy)         (2yz-2sx)     (s2-x2-y2+z2)         0         |
  |         0                 0                 0         (s2+x2+y2+z2) |
  +-                                                                  -+