@ \section{Quaternion Arithmetic in {\tt C++}}

@ We implement non-rigorous quaternion arithmetic in {\tt C++}.

@ \subsection{Programming considerations}

@ \subsubsection{Header and Implementation files}

@ Header file.
<<header>>=
#ifndef QUAT_H
#define QUAT_H
<<global system header files>>
<<class declarations>>
<<inline functions>>
#endif
<<local system header files>>=
#include     <iostream.h>	// stream library
#include     <math.h>		// math routines
#include     <stdlib.h>		// standard lib for exit() function

@ Implementation file.
<<*>>=
<<local system header files>>
#include "quat.h"

@ \subsubsection{Operator Overloading}

We use the overloading mechanism to implement our arithmetic operators.
(Overloading is a way of mimicking polymorphism by defining different
forms of operators and functions and allowing the choice of which to use
to be made at compile-time.)

We implement the assignment operators {\tt =,+=,-=,*=,/=} as public member
functions.
We then implement the binary operations {\tt +,-,*,/} as derived operators
from the above operators. To increase efficiency, the derived operators
are inlined friends.

@ \subsection{Class declaration}
<<class declarations>>=
class quat
{
  <<friend functions>>
private:
  float r;
  float i;
  float j;
  float k;
public:
  <<public members>>
};

@ \subsection{Public Member Functions}

@ \subsubsection{Constructors}

@ Uninitialised constructor.
<<public members>>=
quat() {}                   // uninitialized

@ Constructor from a float precision real.
<<public members>>=
quat(float rr) : r(rr),i(0),j(0),k(0) {}

@ Constructor from a real and the three imaginary components.
<<public members>>=
quat( float, float, float, float );
<<*>>=
quat::quat( float rr, float ii, float jj, float kk )
: r(rr),i(ii),j(jj),k(kk)
{}

@ Copy constructor.
<<public members>>=
quat(const quat& b): r(b.r),i(b.i),j(b.j),k(b.k) {} // copy constructor

@ \subsubsection{Inspectors}

@ Inspectors to return instance variables.
<<public members>>=
float R() const { return r ;}
float I() const { return i ;}
float J() const { return j ;}
float K() const { return k ;}

@ Norm

We define the `norm' of a quaternion to be the square of its modulus,
i.e.
\[
	\hbox{\tt norm}(q) = q\bar{q} = {|q|}^2.
\]
<<public members>>=
float norm() const;
<<*>>=
float quat::norm() const
{
  return r*r+i*i+j*j+k*k;
}

@ Absolute value.
<<public members>>=
float abs() const;

@ We derive the absolute value from the norm by taking the square root.
<<*>>=
float quat::abs() const
{
  return sqrt(norm());
}

@ \subsubsection{Mutators}

@ Assignment
<<public members>>=
quat& operator=(const quat& b) {r=b.r; i=b.i; j=b.j; k=b.k; return *this;}

@ Negation
<<public members>>=
quat& negate();
<<*>>=
quat& quat::negate()
{
  r = -r;
  i = -i;
  j = -j;
  k = -k;
  return *this;
}

@ Conjugation.
<<public members>>=
quat& conjugate();
<<*>>=
quat& quat::conjugate()
{
  i = -i;
  j = -j;
  k = -k;
  return *this;
}

@ Squaring.
<<public members>>=
quat& square();

@ Note that when we square a quaternion, the cross product
of the vector part, $(i,j,k)\times(i,j,k)$, is zero.
<<*>>=
quat& quat::square()
{
  float temp = 2*r;
  r = r*r-(i*i+j*j+k*k);
  i *= temp;
  j *= temp;
  k *= temp;
  return *this;
}

@ Inversion.
<<public members>>=
quat& invert();

@ We invert the quaternion by dividing its conjugate by its norm.
(Note that the norm of the conjugate is the same as the norm of the number
itself, so that the final line of the routine is correct.)
<<*>>=
quat& quat::invert()
{
  float temp = norm();
  r /= temp;
  i /= -temp;
  j /= -temp;
  k /= -temp;
  return (*this);	// Okay, same norm.
}

@ \subsubsection{Assigment Arithmetic Operations}

@ The usual (binary) arithmetic
operators are derived from the corresponding assignment operations below.

@ Addition Assignment.
<<public members>>=
quat& operator+=(const quat& b);
<<*>>=
quat& quat::operator+=(const quat& b)
{
  r+=b.r; i+=b.i; j+=b.j; k+=b.k; return *this;
}

@ Subtraction Assignment.
<<public members>>=
quat& operator-=(const quat& b);
<<*>>=
quat& quat::operator-=(const quat& b)
{
  r-=b.r; i-=b.i; j-=b.j; k-=b.k; return *this;
}

@ Scalar Multiplication Assignment.
<<public members>>=
quat& operator*=(const float b);
<<*>>=
quat& quat::operator*=(const float b)
{
  r*=b; i*=b; j*=b; k*=b; return *this;
}

@ Scalar Division Assignment.
<<public members>>=
quat& operator/=(const float b);
@ NB: the scalar had better not be zero!
<<*>>=
quat& quat::operator/=(const float b)
{
  r/=b; i/=b; j/=b; k/=b; return *this;
}

@ Quaternion Multiplication Assignment.

We use a fast algorithm for multiplication. Our method corrects an error
in the algorithm given in
{\tt http://www.dtek.chalmers.se/Datorsys/Project/}
{\tt qjulia/theory/quaternion.html}.
The algorithm should read as follows:

Let $q=(a,b,c,d)$ and $p=(x,y,z,w)$. Then we define
\begin{eqnarray*}
t_0 &=& (d-c)(z-w)\\
t_1 &=& (a+b)(x+y)\\
t_2 &=& (a-b)(z+w)\\
t_3 &=& (d+c)(x-y)\\
t_4 &=& (d-b)(y-z)\\
t_5 &=& (d+b)(y+z)\\
t_6 &=& (a+c)(x-w)\\
t_7 &=& (a-c)(x+w)\\
t_8 &=& t_5+t_6+t_7\\
t_9 &=& (t_4+t_8)/2
\end{eqnarray*}
and we claim that the quaternion product $qp$ is given by
\[
qp =
 \left(
  \begin{array}[c]{lcr}
   t_0 &-t_5 &+t_9\\
   t_1 &-t_8 &+t_9\\
   t_2 &-t_7 &+t_9\\
   t_3 &-t_6 &+t_9
  \end{array}
 \right)
\]
as is easily (if tediously) verified.

Notice that this method uses only 8 multiplications (compared to the
16 multiplications used in the `naive' method for quaternion multiplication),
at the expense of using 10 extra additions.

<<public members>>=
quat& operator*=(const quat& b);
<<*>>=
quat& quat::operator*=(const quat& b)
{
  float t0 = (k-j)*(b.j-b.k);
  float t1 = (r+i)*(b.r+b.i);
  float t2 = (r-i)*(b.j+b.k);
  float t3 = (k+j)*(b.r-b.i);
  float t4 = (k-i)*(b.i-b.j);
  float t5 = (k+i)*(b.i+b.j);
  float t6 = (r+j)*(b.r-b.k);
  float t7 = (r-j)*(b.r+b.k);
  float t8 = t5+t6+t7;
  float t9 = (t4+t8)/2;
  r = t0+t9-t5;
  i = t1+t9-t8;
  j = t2+t9-t7;
  k = t3+t9-t6;
  return *this;
}

@ Division
<<public members>>=
quat& operator/=(const quat& b);
@ We implement {\bf right}-division assignment
by using the inversion operation and the multiplication assignment.
Unfortunately, the last line of this code looks rather messy, it may
actually be able to call the multiplication assignment operator directly,
rather than having to de-reference the current object.
<<*>>=
quat& quat::operator/=(const quat& b)
{
  quat temp = b;
  return (*this)*=(temp.invert());	// Ugh!
}

@ \subsection{Friend Functions}

@ \subsubsection{Comparison operators}

@ Equality and inequality are the only operators relevant here.
<<friend functions>>=
friend int operator==(const quat& a, const quat& b);
<<*>>=
int operator==(const quat& a, const quat& b)
{
  return (a.r==b.r && a.i==b.i && a.j==b.j && a.k==b.k);
}

<<friend functions>>=
friend int operator!=(const quat& a, const quat& b);
@ For consistency, we derive inequality from equality by negation.
<<*>>=
int operator!=(const quat& a, const quat& b)
{
  return !(a==b); // derived operator
}

@ \subsubsection{Unary Operators}

@ Unary Minus.
<<friend functions>>=
friend quat operator-(quat& x);
<<*>>=
quat operator-(quat& x)
{
  return x.negate();
}

@ \subsection{Derived Arithmetic Operators}

@ Here we implement the usual binary arithmetic operations, by deriving them
from their assignment counterparts. This is considered to be good practice
because it avoids duplication of code and encourages consistency of
operation. The derived operators are also inlined to avoid the inefficiency
of making a call to a simple derived operator.

@ Addition.
<<friend functions>>=
friend quat operator+(const quat& x, const quat& y);
<<inline functions>>=
inline quat operator+(const quat& x, const quat& y)
{
  quat z = x;
  return z += y; // derived operator
}

@ Subtraction.
<<friend functions>>=
friend quat operator-(const quat& x, const quat& y);
<<inline functions>>=
inline quat operator-(const quat& x, const quat& y)
{
  quat z = x;
  return z -= y; // derived operator
}

@ (Right) Scalar Multiplication.
<<friend functions>>=
friend quat operator*(const quat& x, const float y);
<<inline functions>>=
inline quat operator*(const quat& x, const float y)
{
  quat z = x;
  return z *= y; // derived operator
}

@ Quaternion Multiplication.
<<friend functions>>=
friend quat operator*(const quat& x, const quat& y);
<<inline functions>>=
inline quat operator*(const quat& x, const quat& y)
{
  quat z = x;
  return z *= y; // derived operator
}

@ Quaternion Division.
<<friend functions>>=
friend quat operator/(const quat& x, const quat& y);
<<inline functions>>=
inline quat operator/(const quat& x, const quat& y)
{
  quat z = x;
  return z /= y; // derived operator
}

@ Square Root - Not yet implemented.
<<friend functions>>=
friend quat sqrt(const quat& x);
<<*>>=
quat sqrt(const quat& x)
{
  cerr<<"sqrt(const quat&) undefined because quat::sqrt() undefined."<<endl;
  exit(1);	// Should be derived from member function.
}

@ \subsection{I/O Stream Functions}

<<global system header files>>=
#include <iostream.h>	// stream library

<<friend functions>>=
friend ostream& operator<<(ostream& ost, const quat& x);
<<*>>=
ostream& operator<<(ostream& ost, const quat& x)
{
  return ost<<"( "<<x.r<<", "<<x.i<<", "<<x.j<<", "<<x.k<<" )";
}