Quaternion type

From Rosetta Code
Revision as of 10:49, 8 October 2024 by RobertE (talk | contribs) (Formatted the link for the Julia implementation a bit better)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Task
Quaternion type
You are encouraged to solve this task according to the task description, using any language you may know.

Quaternions   are an extension of the idea of   complex numbers.

A complex number has a real and complex part,   sometimes written as   a + bi,
where   a   and   b   stand for real numbers, and   i   stands for the square root of minus 1.

An example of a complex number might be   -3 + 2i,  
where the real part,   a   is   -3.0   and the complex part,   b   is   +2.0.

A quaternion has one real part and three imaginary parts,   i,   j,   and   k.

A quaternion might be written as   a + bi + cj + dk.

In the quaternion numbering system:

  •   i∙i = j∙j = k∙k = i∙j∙k = -1,       or more simply,
  •   ii  = jj  = kk  = ijk   = -1.

The order of multiplication is important, as, in general, for two quaternions:

  q1   and   q2:     q1q2 ≠ q2q1.

An example of a quaternion might be   1 +2i +3j +4k

There is a list form of notation where just the numbers are shown and the imaginary multipliers   i,   j,   and   k   are assumed by position.

So the example above would be written as   (1, 2, 3, 4)


Task

Given the three quaternions and their components:

   q  = (1, 2, 3, 4) = (a,  b,  c,  d)
   q1 = (2, 3, 4, 5) = (a1, b1, c1, d1)
   q2 = (3, 4, 5, 6) = (a2, b2, c2, d2) 

And a wholly real number   r = 7.


Create functions   (or classes)   to perform simple maths with quaternions including computing:

  1. The norm of a quaternion:
  2. The negative of a quaternion:
    = (-a, -b, -c, -d)
  3. The conjugate of a quaternion:
    = ( a, -b, -c, -d)
  4. Addition of a real number   r   and a quaternion   q:
    r + q = q + r = (a+r, b, c, d)
  5. Addition of two quaternions:
    q1 + q2 = (a1+a2, b1+b2, c1+c2, d1+d2)
  6. Multiplication of a real number and a quaternion:
    qr = rq = (ar, br, cr, dr)
  7. Multiplication of two quaternions   q1   and   q2   is given by:
    ( a1a2 − b1b2 − c1c2 − d1d2,
      a1b2 + b1a2 + c1d2 − d1c2,
      a1c2 − b1d2 + c1a2 + d1b2,
      a1d2 + b1c2 − c1b2 + d1a2 )
  8. Show that, for the two quaternions   q1   and   q2:
    q1q2 ≠ q2q1


If a language has built-in support for quaternions, then use it.


C.f.
  •   Vector products
  •   On Quaternions;   or on a new System of Imaginaries in Algebra.   By Sir William Rowan Hamilton LL.D, P.R.I.A., F.R.A.S., Hon. M. R. Soc. Ed. and Dub., Hon. or Corr. M. of the Royal or Imperial Academies of St. Petersburgh, Berlin, Turin and Paris, Member of the American Academy of Arts and Sciences, and of other Scientific Societies at Home and Abroad, Andrews' Prof. of Astronomy in the University of Dublin, and Royal Astronomer of Ireland.



Action!

INCLUDE "H6:REALMATH.ACT"

DEFINE A_="+0"
DEFINE B_="+6"
DEFINE C_="+12"
DEFINE D_="+18"

TYPE Quaternion=[CARD a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3]
REAL neg

PROC Init()
  ValR("-1",neg)
RETURN

BYTE FUNC Positive(REAL POINTER x)
  BYTE ARRAY tmp

  tmp=x
  IF (tmp(0)&$80)=$00 THEN
    RETURN (1)
  FI
RETURN (0)

PROC PrintQuat(Quaternion POINTER q)
  PrintR(q A_)
  IF Positive(q B_) THEN Put('+) FI
  PrintR(q B_) Put('i)
  IF Positive(q C_) THEN Put('+) FI
  PrintR(q C_) Put('j)
  IF Positive(q D_) THEN Put('+) FI
  PrintR(q D_) Put('k)
RETURN

PROC PrintQuatE(Quaternion POINTER q)
  PrintQuat(q) PutE()
RETURN

PROC QuatIntInit(Quaternion POINTER q INT ia,ib,ic,id)
  IntToReal(ia,q A_)
  IntToReal(ib,q B_)
  IntToReal(ic,q C_)
  IntToReal(id,q D_)
RETURN

PROC Sqr(REAL POINTER a,b)
  RealMult(a,a,b)
RETURN

PROC QuatNorm(Quaternion POINTER q REAL POINTER res)
  REAL r1,r2,r3

  Sqr(q A_,r1)      ;r1=q.a^2
  Sqr(q B_,r2)      ;r2=q.b^2
  RealAdd(r1,r2,r3) ;r3=q.a^2+q.b^2
  Sqr(q C_,r1)      ;r1=q.c^2
  RealAdd(r3,r1,r2) ;r2=q.a^2+q.b^2+q.c^2
  Sqr(q D_,r1)      ;r1=q.d^2
  RealAdd(r2,r1,r3) ;r3=q.a^2+q.b^2+q.c^2+q.d^2
  Sqrt(r3,res)      ;res=sqrt(q.a^2+q.b^2+q.c^2+q.d^2)
RETURN

PROC QuatNegative(Quaternion POINTER q,res)
  RealMult(q A_,neg,res A_) ;res.a=-q.a
  RealMult(q B_,neg,res B_) ;res.b=-q.b
  RealMult(q C_,neg,res C_) ;res.c=-q.c
  RealMult(q D_,neg,res D_) ;res.d=-q.d
RETURN

PROC QuatConjugate(Quaternion POINTER q,res)
  RealAssign(q A_,res A_)   ;res.a=q.a
  RealMult(q B_,neg,res B_) ;res.b=-q.b
  RealMult(q C_,neg,res C_) ;res.c=-q.c
  RealMult(q D_,neg,res D_) ;res.d=-q.d
RETURN

PROC QuatAddReal(Quaternion POINTER q REAL POINTER r
                 Quaternion POINTER res)
  RealAdd(q A_,r,res A_)  ;res.a=q.a+r
  RealAssign(q B_,res B_) ;res.b=q.b
  RealAssign(q C_,res C_) ;res.c=q.c
  RealAssign(q D_,res D_) ;res.d=q.d
RETURN

PROC QuatAdd(Quaternion POINTER q1,q2,res)
  RealAdd(q1 A_,q2 A_,res A_) ;res.a=q1.a+q2.a
  RealAdd(q1 B_,q2 B_,res B_) ;res.b=q1.b+q2.b
  RealAdd(q1 C_,q2 C_,res C_) ;res.c=q1.c+q2.c
  RealAdd(q1 D_,q2 D_,res D_) ;res.d=q1.d+q2.d
RETURN

PROC QuatMultReal(Quaternion POINTER q REAL POINTER r
                  Quaternion POINTER res)
  RealMult(q A_,r,res A_) ;res.a=q.a*r
  RealMult(q B_,r,res B_) ;res.b=q.b*r
  RealMult(q C_,r,res C_) ;res.c=q.c*r
  RealMult(q D_,r,res D_) ;res.d=q.d*r
RETURN

PROC QuatMult(Quaternion POINTER q1,q2,res)
  REAL r1,r2

  RealMult(q1 A_,q2 A_,r1) ;r1=q1.a*q2.a
  RealMult(q1 B_,q2 B_,r2) ;r2=q1.b*q2.b
  RealSub(r1,r2,r3)        ;r3=q1.a*q2.a-q1.b*q2.b
  RealMult(q1 C_,q2 C_,r1) ;r1=q1.c*q2.c
  RealSub(r3,r1,r2)        ;r2=q1.a*q2.a-q1.b*q2.b-q1.c*q2.c
  RealMult(q1 D_,q2 D_,r1) ;r1=q1.d*q2.d
  RealSub(r2,r1,res A_)    ;res.a=q1.a*q2.a-q1.b*q2.b-q1.c*q2.c-q1.d*q2.d

  RealMult(q1 A_,q2 B_,r1) ;r1=q1.a*q2.b
  RealMult(q1 B_,q2 A_,r2) ;r2=q1.b*q2.a
  RealAdd(r1,r2,r3)        ;r3=q1.a*q2.b+q1.b*q2.a
  RealMult(q1 C_,q2 D_,r1) ;r1=q1.c*q2.d
  RealAdd(r3,r1,r2)        ;r2=q1.a*q2.b+q1.b*q2.a+q1.c*q2.d
  RealMult(q1 D_,q2 C_,r1) ;r1=q1.d*q2.c
  RealSub(r2,r1,res B_)    ;res.b=q1.a*q2.b+q1.b*q2.a+q1.c*q2.d-q1.d*q2.c

  RealMult(q1 A_,q2 C_,r1) ;r1=q1.a*q2.c
  RealMult(q1 B_,q2 D_,r2) ;r2=q1.b*q2.d
  RealSub(r1,r2,r3)        ;r3=q1.a*q2.c-q1.b*q2.d
  RealMult(q1 C_,q2 A_,r1) ;r1=q1.c*q2.a
  RealAdd(r3,r1,r2)        ;r2=q1.a*q2.c-q1.b*q2.d+q1.c*q2.a
  RealMult(q1 D_,q2 B_,r1) ;r1=q1.d*q2.b
  RealAdd(r2,r1,res C_)    ;res.c=q1.a*q2.c-q1.b*q2.d+q1.c*q2.a+q1.d*q2.b

  RealMult(q1 A_,q2 D_,r1) ;r1=q1.a*q2.d
  RealMult(q1 B_,q2 C_,r2) ;r2=q1.b*q2.c
  RealAdd(r1,r2,r3)        ;r3=q1.a*q2.d+q1.b*q2.c
  RealMult(q1 C_,q2 B_,r1) ;r1=q1.c*q2.b
  RealSub(r3,r1,r2)        ;r2=q1.a*q2.d+q1.b*q2.c-q1.c*q2.b
  RealMult(q1 D_,q2 A_,r1) ;r1=q1.d*q2.a
  RealAdd(r2,r1,res D_)    ;res.d=q1.a*q2.d+q1.b*q2.c-q1.c*q2.b+q1.d*q2.a
RETURN

PROC Main()
  Quaternion q,q1,q2,q3
  REAL r,r2

  Put(125) PutE() ;clear the screen
  MathInit()
  Init()

  QuatIntInit(q,1,2,3,4)
  QuatIntInit(q1,2,3,4,5)
  QuatIntInit(q2,3,4,5,6)
  IntToReal(7,r)

  Print(" q = ") PrintQuatE(q)
  Print("q1 = ") PrintQuatE(q1)
  Print("q2 = ") PrintQuatE(q2)
  Print(" r = ") PrintRE(r) PutE()

  QuatNorm(q,r2) Print(" Norm(q) = ") PrintRE(r2)
  QuatNorm(q1,r2) Print("Norm(q1) = ") PrintRE(r2)
  QuatNorm(q2,r2) Print("Norm(q2) = ") PrintRE(r2)
  QuatNegative(q,q3) Print("      -q = ") PrintQuatE(q3)
  QuatConjugate(q,q3) Print(" Conj(q) = ") PrintQuatE(q3)
  QuatAddReal(q,r,q3) Print("     q+r = ") PrintQuatE(q3)
  QuatAdd(q1,q2,q3) Print("   q1+q2 = ") PrintQuatE(q3)
  QuatAdd(q2,q1,q3) Print("   q2+q1 = ") PrintQuatE(q3)
  QuatMultReal(q,r,q3) Print("     q*r = ") PrintQuatE(q3)
  QuatMult(q1,q2,q3) Print("   q1*q2 = ") PrintQuatE(q3)
  QuatMult(q2,q1,q3) Print("   q2*q1 = ") PrintQuatE(q3)
RETURN
Output:

Screenshot from Atari 8-bit computer

 q = 1+2i+3j+4k
q1 = 2+3i+4j+5k
q2 = 3+4i+5j+6k
 r = 7

 Norm(q) = 5.47722543
Norm(q1) = 7.34846906
Norm(q2) = 9.27361833
      -q = -1-2i-3j-4k
 Conj(q) = 1-2i-3j-4k
     q+r = 8+2i+3j+4k
   q1+q2 = 5+7i+9j+11k
   q2+q1 = 5+7i+9j+11k
     q*r = 7+14i+21j+28k
   q1*q2 = -56+16i+24j+26k
   q2*q1 = -56+18i+20j+28k

Ada

The package specification (works with any floating-point type):

generic
   type Real is digits <>;
package Quaternions is
   type Quaternion is record
      A, B, C, D : Real;
   end record;
   function "abs" (Left : Quaternion) return Real;
   function Conj (Left : Quaternion) return Quaternion;
   function "-" (Left : Quaternion) return Quaternion;
   function "+" (Left, Right : Quaternion) return Quaternion;
   function "-" (Left, Right : Quaternion) return Quaternion;
   function "*" (Left : Quaternion; Right : Real) return Quaternion;
   function "*" (Left : Real; Right : Quaternion) return Quaternion;
   function "*" (Left, Right : Quaternion) return Quaternion;
   function Image (Left : Quaternion) return String;
end Quaternions;

The package implementation:

with Ada.Numerics.Generic_Elementary_Functions;
package body Quaternions is
   package Elementary_Functions is
      new Ada.Numerics.Generic_Elementary_Functions (Real);
   use Elementary_Functions;
   function "abs" (Left : Quaternion) return Real is
   begin
      return Sqrt (Left.A**2 + Left.B**2 + Left.C**2 + Left.D**2);
   end "abs";
   function Conj (Left : Quaternion) return Quaternion is
   begin
      return (A => Left.A, B => -Left.B, C => -Left.C, D => -Left.D);
   end Conj;
   function "-" (Left : Quaternion) return Quaternion is
   begin
      return (A => -Left.A, B => -Left.B, C => -Left.C, D => -Left.D);
   end "-";
   function "+" (Left, Right : Quaternion) return Quaternion is
   begin
      return
      (  A => Left.A + Right.A, B => Left.B + Right.B, 
         C => Left.C + Right.C, D => Left.D + Right.D
      );
   end "+";
   function "-" (Left, Right : Quaternion) return Quaternion is
   begin
      return
      (  A => Left.A - Right.A, B => Left.B - Right.B, 
         C => Left.C - Right.C, D => Left.D - Right.D
      );
   end "-";
   function "*" (Left : Quaternion; Right : Real) return Quaternion is
   begin
      return
      (  A => Left.A * Right, B => Left.B * Right, 
         C => Left.C * Right, D => Left.D * Right
      );
   end "*";
   function "*" (Left : Real; Right : Quaternion) return Quaternion is
   begin
      return Right * Left;
   end "*";
   function "*" (Left, Right : Quaternion) return Quaternion is
   begin
      return
      (  A => Left.A * Right.A - Left.B * Right.B - Left.C * Right.C - Left.D * Right.D,
         B => Left.A * Right.B + Left.B * Right.A + Left.C * Right.D - Left.D * Right.C,
         C => Left.A * Right.C - Left.B * Right.D + Left.C * Right.A + Left.D * Right.B,
         D => Left.A * Right.D + Left.B * Right.C - Left.C * Right.B + Left.D * Right.A
      );
   end "*";
   function Image (Left : Quaternion) return String is
   begin
      return Real'Image (Left.A) & " +"  &
             Real'Image (Left.B) & "i +" &
             Real'Image (Left.C) & "j +" &
             Real'Image (Left.D) & "k";
   end Image;
end Quaternions;

Test program:

with Ada.Text_IO;  use Ada.Text_IO;
with Quaternions;
procedure Test_Quaternion is
   package Float_Quaternion is new Quaternions (Float);
   use Float_Quaternion;
   q  : Quaternion := (1.0, 2.0, 3.0, 4.0);
   q1 : Quaternion := (2.0, 3.0, 4.0, 5.0);
   q2 : Quaternion := (3.0, 4.0, 5.0, 6.0);
   r  : Float      := 7.0;
begin
   Put_Line ("q = "       & Image (q));
   Put_Line ("q1 = "      & Image (q1));
   Put_Line ("q2 = "      & Image (q2));
   Put_Line ("r ="        & Float'Image (r));
   Put_Line ("abs q ="    & Float'Image (abs q));
   Put_Line ("abs q1 ="   & Float' Image (abs q1));
   Put_Line ("abs q2 ="   & Float' Image (abs q2));
   Put_Line ("-q = "      & Image (-q));
   Put_Line ("conj q = "  & Image (Conj (q)));
   Put_Line ("q1 + q2 = " & Image (q1 + q2));
   Put_Line ("q2 + q1 = " & Image (q2 + q1));
   Put_Line ("q * r = "   & Image (q * r));
   Put_Line ("r * q = "   & Image (r * q));
   Put_Line ("q1 * q2 = " & Image (q1 * q2));
   Put_Line ("q2 * q1 = " & Image (q2 * q1));
end Test_Quaternion;
Output:
q =  1.00000E+00 + 2.00000E+00i + 3.00000E+00j + 4.00000E+00k
q1 =  2.00000E+00 + 3.00000E+00i + 4.00000E+00j + 5.00000E+00k
q2 =  3.00000E+00 + 4.00000E+00i + 5.00000E+00j + 6.00000E+00k
r = 7.00000E+00
abs q = 5.47723E+00
abs q1 = 7.34847E+00
abs q2 = 9.27362E+00
-q = -1.00000E+00 +-2.00000E+00i +-3.00000E+00j +-4.00000E+00k
conj q =  1.00000E+00 +-2.00000E+00i +-3.00000E+00j +-4.00000E+00k
q1 + q2 =  5.00000E+00 + 7.00000E+00i + 9.00000E+00j + 1.10000E+01k
q2 + q1 =  5.00000E+00 + 7.00000E+00i + 9.00000E+00j + 1.10000E+01k
q * r =  7.00000E+00 + 1.40000E+01i + 2.10000E+01j + 2.80000E+01k
r * q =  7.00000E+00 + 1.40000E+01i + 2.10000E+01j + 2.80000E+01k
q1 * q2 = -5.60000E+01 + 1.60000E+01i + 2.40000E+01j + 2.60000E+01k
q2 * q1 = -5.60000E+01 + 1.80000E+01i + 2.00000E+01j + 2.80000E+01k

ALGOL 68

Translation of: python

Note: This specimen retains the original python coding style.

Works with: ALGOL 68 version Revision 1 - one minor extension to language used - PRAGMA READ, similar to C's #include directive.
Works with: ALGOL 68G version Any - tested with release algol68g-2.6.

File: prelude/Quaternion.a68

# -*- coding: utf-8 -*- #

COMMENT REQUIRES:
  MODE QUATSCAL = REAL; # Scalar #
  QUATSCAL quat small scal = small real;
END COMMENT

# PROVIDES: #
FORMAT quat scal fmt := $g(-0, 4)$;
FORMAT signed fmt = $b("+", "")f(quat scal fmt)$;

FORMAT quat fmt = $f(quat scal fmt)"+"f(quat scal fmt)"i+"f(quat scal fmt)"j+"f(quat scal fmt)"k"$;
FORMAT squat fmt = $f(signed fmt)f(signed fmt)"i"f(signed fmt)"j"f(signed fmt)"k"$;

MODE QUAT = STRUCT(QUATSCAL r, i, j, k);
QUAT i=(0, 1, 0, 0),
     j=(0, 0, 1, 0),
     k=(0, 0, 0, 1);

MODE QUATCOSCAL = UNION(INT, SHORT REAL, SHORT INT); 
MODE QUATSUBSCAL = UNION(QUATCOSCAL, QUATSCAL);

MODE COMPLSCAL = STRUCT(QUATSCAL r, im);
# compatable but not the same #
MODE ISOQUAT = UNION([]REAL, []INT, []SHORT REAL, []SHORT INT, []QUATSCAL);
MODE COQUAT  = UNION(COMPLSCAL, QUATCOSCAL, ISOQUAT);
MODE SUBQUAT = UNION(COQUAT, QUAT); # subset is itself #

MODE QUATERNION = QUAT;

PROC quat fix type error = (QUAT quat, []STRING msg)BOOL: (
  putf(stand error, ($"Type error:"$,$" "g$, msg, quat fmt, quat, $l$));
  stop
);

COMMENT
For a list of coercions expected in A68 c.f. 
* http://rosettacode.org/wiki/ALGOL_68#Coercion_.28casting.29 # ...

Pre-Strong context: Deproceduring, dereferencing & uniting. e.g. OP arguments
  * soft(deproceduring for assignment), 
  * weak(dereferencing for slicing and OF selection), 
  * meek(dereferencing for indexing, enquiries and PROC calls),
  * firm(uniting of OPerators), 
Strong context only: widening (INT=>REAL=>COMPL), rowing (REAL=>[]REAL) & voiding
  * strong(widening,rowing,voiding for identities/initialisations, arguments and casts et al)
Key points:
  * arguments to OPerators do not widen or row!
  * UNITING is permitted in OP/String ccontext.

There are 4 principle scenerios for most operators:
+---------------+-------------------------------+-------------------------------+
|  OP e.g. *    |  SCALar                       |  QUATernion                   |
+---------------+-------------------------------+-------------------------------+
|  SCALar       |  SCAL * SCAL ... inherit      |  SCAL * QUAT                  |
+---------------+-------------------------------+-------------------------------+
|  QUATernion   |  QUAT * SCAL                  |  QUAT * QUAT                  |
+---------------+-------------------------------+-------------------------------+
However this is compounded with SUBtypes of the SCALar & isomorphs the QUATernion, 
e.g. 
* SCAL may be a superset of SHORT REAL or INT - a widening coercion is required
* QUAT may be a superset eg of COMPL or [4]INT
* QUAT may be a structural isomorph eg of [4]REAL
+---------------+---------------+---------------+---------------+---------------+
|  OP e.g. *    |  SUBSCAL      |  SCALar       |  COQUAT       |  QUATernion   |
+---------------+---------------+---------------+---------------+---------------+
|  SUBSCAL      |                               |  inherit      |  SUBSCAT*QUAT |
+---------------+           inherit             +---------------+---------------+
|  SCALar       |                               |  inherit      |  SCAL * QUAT  |
+---------------+---------------+---------------+---------------+---------------+
|  COQUAT       |  inherit      |  inherit      |  inherit      |  COQUAT*QUAT  |
+---------------+---------------+---------------+---------------+---------------+
|  QUATernion   | QUAT*SUBSCAL  |  QUAT*SCAL    | QUAT * COQUAT |  QUAT * QUAT  |
+---------------+---------------+---------------+---------------+---------------+
Keypoint: if an EXPLICIT QUAT is not involved, then we can simple inherit, OR QUATINIT!
END COMMENT

MODE CLASSQUAT = STRUCT(
    PROC (REF QUAT #new#, QUATSCAL #r#, QUATSCAL #i#, QUATSCAL #j#, QUATSCAL #k#)REF QUAT new,
    PROC (REF QUAT #self#)QUAT conjugate,
    PROC (REF QUAT #self#)QUATSCAL norm sq,
    PROC (REF QUAT #self#)QUATSCAL norm,
    PROC (REF QUAT #self#)QUAT reciprocal,
    PROC (REF QUAT #self#)STRING repr,
    PROC (REF QUAT #self#)QUAT neg,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT add,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT radd,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT sub,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT mul,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT rmul,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT div,
    PROC (REF QUAT #self#, SUBQUAT #other#)QUAT rdiv,
    PROC (REF QUAT #self#)QUAT exp
);

CLASSQUAT class quat = (

  # PROC new =#(REF QUAT new, QUATSCAL r, i, j, k)REF QUAT: (
        # 'Defaults all parts of quaternion to zero' #
        IF new ISNT REF QUAT(NIL) THEN new ELSE HEAP QUAT FI := (r, i, j, k)
    ),

  # PROC conjugate =#(REF QUAT self)QUAT:
        (r OF self, -i OF self, -j OF self, -k OF self),

  # PROC norm sq =#(REF QUAT self)QUATSCAL:
        r OF self**2 + i OF self**2 + j OF self**2 + k OF self**2,

  # PROC norm =#(REF QUAT self)QUATSCAL:
        sqrt((norm sq OF class quat)(self)),

  # PROC reciprocal =#(REF QUAT self)QUAT:(
        QUATSCAL n2 = (norm sq OF class quat)(self);
        QUAT conj = (conjugate OF class quat)(self);
        (r OF conj/n2, i OF conj/n2, j OF conj/n2, k OF conj/n2)
    ),

  # PROC repr =#(REF QUAT self)STRING: (
        # 'Shorter form of Quaternion as string' #
        FILE f; STRING s; associate(f, s);
        putf(f, (squat fmt, r OF self>=0, r OF self,
             i OF self>=0, i OF self, j OF self>=0, j OF self, k OF self>=0, k OF self));
        close(f);
        s
    ),

  # PROC neg =#(REF QUAT self)QUAT:
        (-r OF self, -i OF self, -j OF self, -k OF self),

  # PROC add =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
            (QUAT other): (r OF self + r OF other, i OF self + i OF other, j OF self + j OF other, k OF self + k OF other),
            (QUATSUBSCAL other): (r OF self + QUATSCALINIT other, i OF self, j OF self, k OF self)
        OUT IF quat fix type error(SKIP,"in add") THEN SKIP ELSE stop FI
        ESAC,

  # PROC radd =#(REF QUAT self, SUBQUAT other)QUAT:
        (add OF class quat)(self, other),

  # PROC sub =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
            (QUAT other): (r OF self - r OF other, i OF self - i OF other, j OF self - j OF other, k OF self - k OF other),
            (QUATSCAL other): (r OF self - other, i OF self, j OF self, k OF self)
        OUT IF quat fix type error(self,"in sub") THEN SKIP ELSE stop FI
        ESAC,

  # PROC mul =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
            (QUAT other):(
                 r OF self*r OF other - i OF self*i  OF other - j OF self*j  OF other - k OF self*k  OF other,
                 r OF self*i  OF other + i OF self*r OF other + j OF self*k  OF other - k OF self*j  OF other,
                 r OF self*j  OF other - i OF self*k  OF other + j OF self*r OF other + k OF self*i  OF other,
                 r OF self*k  OF other + i OF self*j  OF other - j OF self*i  OF other + k OF self*r OF other
            ),
            (QUATSCAL other): ( r OF self * other, i OF self * other, j OF self * other, k OF self * other)
        OUT IF quat fix type error(self,"in mul") THEN SKIP ELSE stop FI
        ESAC,

  # PROC rmul =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
          (QUAT other): (mul OF class quat)(LOC QUAT := other, self),
          (QUATSCAL other): (mul OF class quat)(self, other)
        OUT IF quat fix type error(self,"in rmul") THEN SKIP ELSE stop FI
        ESAC,

  # PROC div =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
            (QUAT other): (mul OF class quat)(self, (reciprocal OF class quat)(LOC QUAT := other)),
            (QUATSCAL other): (mul OF class quat)(self, 1/other)
        OUT IF quat fix type error(self,"in div") THEN SKIP ELSE stop FI
        ESAC,

  # PROC rdiv =#(REF QUAT self, SUBQUAT other)QUAT:
        CASE other IN
          (QUAT other): (div OF class quat)(LOC QUAT := other, self),
          (QUATSCAL other): (div OF class quat)(LOC QUAT := (other, 0, 0, 0), self)
        OUT IF quat fix type error(self,"in rdiv") THEN SKIP ELSE stop FI
        ESAC,

  # PROC exp =#(REF QUAT self)QUAT: (
    QUAT fac := self;
    QUAT sum := 1.0 + fac;
    FOR i FROM 2 TO bits width WHILE ABS(fac + quat small scal) /= quat small scal DO
      VOID(sum +:= (fac *:= self / ##QUATSCAL(i)))
    OD;
    sum
  )
);

PRIO INIT = 1;
OP QUATSCALINIT = (QUATSUBSCAL scal)QUATSCAL: 
  CASE scal IN
    (INT scal): scal,
    (SHORT INT scal): scal,
    (SHORT REAL scal): scal
    OUT IF quat fix type error(SKIP,"in QUATSCALINIT") THEN SKIP ELSE stop FI
  ESAC;

OP INIT = (REF QUAT new, SUBQUAT from)REF QUAT:
  new := 
    CASE from IN
      (QUATSUBSCAL scal):(QUATSCALINIT scal, 0, 0, 0)
      #(COQUAT rijk):(new OF class quat)(LOC QUAT := new, rijk[1], rijk[2], rijk[3], rijk[4]),#
    OUT IF quat fix type error(SKIP,"in INIT") THEN SKIP ELSE stop FI
    ESAC;


OP QUATINIT = (COQUAT lhs)REF QUAT: (HEAP QUAT)INIT lhs;

OP +    = (QUAT q)QUAT:   q,
   -    = (QUAT q)QUAT:   (neg  OF class quat)(LOC QUAT := q),
   CONJ = (QUAT q)QUAT:   (conjugate OF class quat)(LOC QUAT := q),
   ABS  = (QUAT q)QUATSCAL:   (norm OF class quat)(LOC QUAT := q),
   REPR = (QUAT q)STRING: (repr OF class quat)(LOC QUAT := q);
# missing: Diadic: I, J, K END #

OP +:= = (REF QUAT a, QUAT b)QUAT: a:=( add OF class quat)(a, b),
   +:= = (REF QUAT a, COQUAT b)QUAT: a:=( add OF class quat)(a, b),
   +=: = (QUAT a, REF QUAT b)QUAT: b:=(radd OF class quat)(b, a),
   +=: = (COQUAT a, REF QUAT b)QUAT: b:=(radd OF class quat)(b, a);
# missing: Worthy PLUSAB, PLUSTO for SHORT/LONG INT QUATSCAL & COMPL #

OP -:= = (REF QUAT a, QUAT b)QUAT: a:=( sub OF class quat)(a, b),
   -:= = (REF QUAT a, COQUAT b)QUAT: a:=( sub OF class quat)(a, b);
# missing: Worthy MINUSAB for SHORT/LONG INT ##COQUAT & COMPL #

PRIO *=: = 1, /=: = 1;
OP *:= = (REF QUAT a, QUAT b)QUAT: a:=( mul OF class quat)(a, b),
   *:= = (REF QUAT a, COQUAT b)QUAT: a:=( mul OF class quat)(a, b),
   *=: = (QUAT a, REF QUAT b)QUAT: b:=(rmul OF class quat)(b, a),
   *=: = (COQUAT a, REF QUAT b)QUAT: b:=(rmul OF class quat)(b, a);
# missing: Worthy TIMESAB, TIMESTO for SHORT/LONG INT ##COQUAT & COMPL #

OP /:= = (REF QUAT a, QUAT b)QUAT: a:=( div OF class quat)(a, b),
   /:= = (REF QUAT a, COQUAT b)QUAT: a:=( div OF class quat)(a, b),
   /=: = (QUAT a, REF QUAT b)QUAT: b:=(rdiv OF class quat)(b, a),
   /=: = (COQUAT a, REF QUAT b)QUAT: b:=(rdiv OF class quat)(b, a);
# missing: Worthy OVERAB, OVERTO for SHORT/LONG INT ##COQUAT & COMPL #

OP + = (QUAT a, b)QUAT:      ( add OF class quat)(LOC QUAT := a, b),
   + = (QUAT a, COQUAT b)QUAT: ( add OF class quat)(LOC QUAT := a, b),
   + = (COQUAT a, QUAT b)QUAT: (radd OF class quat)(LOC QUAT := b, a);
 
OP - = (QUAT a, b)QUAT:      ( sub OF class quat)(LOC QUAT := a, b),
   - = (QUAT a, COQUAT b)QUAT: ( sub OF class quat)(LOC QUAT := a, b),
   - = (COQUAT a, QUAT b)QUAT:-( sub OF class quat)(LOC QUAT := b, a);
 
OP * = (QUAT a, b)QUAT:      ( mul OF class quat)(LOC QUAT := a, b),
   * = (QUAT a, COQUAT b)QUAT: ( mul OF class quat)(LOC QUAT := a, b),
   * = (COQUAT a, QUAT b)QUAT: (rmul OF class quat)(LOC QUAT := b, a);
 
OP / = (QUAT a, b)QUAT:      ( div OF class quat)(LOC QUAT := a, b),
   / = (QUAT a, COQUAT b)QUAT: ( div OF class quat)(LOC QUAT := a, b),
   / = (COQUAT a, QUAT b)QUAT: 
         ( div OF class quat)(LOC QUAT := QUATINIT 1, a);

PROC quat exp = (QUAT q)QUAT:   (exp OF class quat)(LOC QUAT := q);

SKIP # missing: quat arc{sin, cos, tan}h, log, exp, ln etc END #

File: test/Quaternion.a68

#!/usr/bin/a68g --script #
# -*- coding: utf-8 -*- #

# REQUIRES: #
  MODE QUATSCAL = REAL; # Scalar #
  QUATSCAL quat small scal = small real;

PR READ "prelude/Quaternion.a68" PR;

test:(
    REAL r = 7;
    QUAT q  = (1, 2, 3, 4),
         q1 = (2, 3, 4, 5),
         q2 = (3, 4, 5, 6);

    printf((
        $"r = "      f(quat scal fmt)l$, r,
        $"q = "      f(quat fmt)l$, q,
        $"q1 = "     f(quat fmt)l$, q1,
        $"q2 = "     f(quat fmt)l$, q2,
        $"ABS q = "  f(quat scal fmt)", "$, ABS q,
        $"ABS q1 = " f(quat scal fmt)", "$, ABS q1,
        $"ABS q2 = " f(quat scal fmt)l$, ABS q2,
        $"-q = "     f(quat fmt)l$, -q,
        $"CONJ q = " f(quat fmt)l$, CONJ q,
        $"r + q = "  f(quat fmt)l$, r + q,
        $"q + r = "  f(quat fmt)l$, q + r,
        $"q1 + q2 = "f(quat fmt)l$, q1 + q2,
        $"q2 + q1 = "f(quat fmt)l$, q2 + q1,
        $"q * r = "  f(quat fmt)l$, q * r,
        $"r * q = "  f(quat fmt)l$, r * q,
        $"q1 * q2 = "f(quat fmt)l$, q1 * q2,
        $"q2 * q1 = "f(quat fmt)l$, q2 * q1
    ));

CO
        $"ASSERT q1 * q2 != q2 * q1 = "f(quat fmt)l$, ASSERT q1 * q2 != q2 * q1, $l$;
END CO

    printf((
        $"i*i = "         f(quat fmt)l$, i*i,
        $"j*j = "         f(quat fmt)l$, j*j,
        $"k*k = "         f(quat fmt)l$, k*k,
        $"i*j*k = "       f(quat fmt)l$, i*j*k,
        $"q1 / q2 = "     f(quat fmt)l$, q1 / q2,
        $"q1 / q2 * q2 = "f(quat fmt)l$, q1 / q2 * q2,
        $"q2 * q1 / q2 = "f(quat fmt)l$, q2 * q1 / q2,
        $"1/q1 * q1 = "   f(quat fmt)l$, 1.0/q1 * q1,
        $"q1 / q1 = "     f(quat fmt)l$, q1 / q1,
        $"quat exp(pi * i) = " f(quat fmt)l$, quat exp(pi * i),
        $"quat exp(pi * j) = " f(quat fmt)l$, quat exp(pi * j),
        $"quat exp(pi * k) = " f(quat fmt)l$, quat exp(pi * k)
    ));
    print((REPR(-q1*q2), ", ", REPR(-q2*q1), new line))
)
Output:
r = 7.0000
q = 1.0000+2.0000i+3.0000j+4.0000k
q1 = 2.0000+3.0000i+4.0000j+5.0000k
q2 = 3.0000+4.0000i+5.0000j+6.0000k
ABS q = 5.4772, ABS q1 = 7.3485, ABS q2 = 9.2736
-q = -1.0000+-2.0000i+-3.0000j+-4.0000k
CONJ q = 1.0000+-2.0000i+-3.0000j+-4.0000k
r + q = 8.0000+2.0000i+3.0000j+4.0000k
q + r = 8.0000+2.0000i+3.0000j+4.0000k
q1 + q2 = 5.0000+7.0000i+9.0000j+11.0000k
q2 + q1 = 5.0000+7.0000i+9.0000j+11.0000k
q * r = 7.0000+14.0000i+21.0000j+28.0000k
r * q = 7.0000+14.0000i+21.0000j+28.0000k
q1 * q2 = -56.0000+16.0000i+24.0000j+26.0000k
q2 * q1 = -56.0000+18.0000i+20.0000j+28.0000k
i*i = -1.0000+.0000i+.0000j+.0000k
j*j = -1.0000+.0000i+.0000j+.0000k
k*k = -1.0000+.0000i+.0000j+.0000k
i*j*k = -1.0000+.0000i+.0000j+.0000k
q1 / q2 = .7907+.0233i+-.0000j+.0465k
q1 / q2 * q2 = 2.0000+3.0000i+4.0000j+5.0000k
q2 * q1 / q2 = 2.0000+3.4651i+3.9070j+4.7674k
1/q1 * q1 = 2.0000+3.0000i+4.0000j+5.0000k
q1 / q1 = 1.0000+.0000i+.0000j+.0000k
quat exp(pi * i) = -1.0000+.0000i+.0000j+.0000k
quat exp(pi * j) = -1.0000+.0000i+.0000j+.0000k
quat exp(pi * k) = -1.0000+.0000i+.0000j+.0000k
+56.0000-16.0000i-24.0000j-26.0000k, +56.0000-18.0000i-20.0000j-28.0000k

ALGOL W

begin
    % Quaternion record type                                                 %
    record Quaternion ( real a, b, c, d );

    % returns the norm of the specified quaternion                           %
    real procedure normQ ( reference(Quaternion) value q ) ;
        sqrt( (a(q) * a(q)) + (b(q) * b(q)) + (c(q) * c(q)) + (d(q) * d(q)) ); 

    % returns the negative of the specified quaternion                       %
    reference(Quaternion) procedure negQ ( reference(Quaternion) value q ) ;
        Quaternion( - a(q), - b(q), - c(q), - d(q) );

    % returns the conjugate of the specified quaternion                      %
    reference(Quaternion) procedure conjQ ( reference(Quaternion) value q ) ;
        Quaternion(   a(q), - b(q), - c(q), - d(q) );

    % returns the sum of a real and a quaternion                             %
    reference(Quaternion) procedure addRQ ( real                  value r
                                          ; reference(Quaternion) value q
                                          ) ;
        Quaternion( r + a(q), b(q), c(q), d(q) );

    % returns the sum of a quaternion and a real                             %
    reference(Quaternion) procedure addQR ( reference(Quaternion) value q
                                          ; real                  value r
                                          ) ;
        Quaternion( r + a(q), b(q), c(q), d(q) );

    % returns the sum of the specified quaternions                           %
    reference(Quaternion) procedure addQQ ( reference(Quaternion) value q1
                                          ; reference(Quaternion) value q2
                                          ) ;
        Quaternion( a(q1) + a(q2), b(q1) + b(q2), c(q1) + c(q2), d(q1) + d(q2) );

    % returns the specified quaternion multiplied by a real                  %
    reference(Quaternion) procedure mulQR ( reference(Quaternion) value q
                                          ; real                  value r
                                          ) ;
        Quaternion( r * a(q), r * b(q), r * c(q), r * d(q) );

    % returns a real multiplied by the specified quaternion                  %
    reference(Quaternion) procedure mulRQ ( real                  value r
                                          ; reference(Quaternion) value q
                                          ) ;
        mulQR( q, r );

    % returns the Quaternion product of the specified quaternions            %
    reference(Quaternion) procedure mulQQ( reference(Quaternion) value q1
                                         ; reference(Quaternion) value q2
                                         ) ;
        Quaternion( (a(q1) * a(q2)) - (b(q1) * b(q2)) - (c(q1) * c(q2)) - (d(q1) * d(q2))
                  , (a(q1) * b(q2)) + (b(q1) * a(q2)) + (c(q1) * d(q2)) - (d(q1) * c(q2))
                  , (a(q1) * c(q2)) - (b(q1) * d(q2)) + (c(q1) * a(q2)) + (d(q1) * b(q2))
                  , (a(q1) * d(q2)) + (b(q1) * c(q2)) - (c(q1) * b(q2)) + (d(q1) * a(q2))
                  );

    % returns true if the two quaternions are equal, false otherwise         %
    logical procedure equalQ( reference(Quaternion) value q1
                            ; reference(Quaternion) value q2
                            ) ;
        a(q1) = a(q2) and b(q1) = b(q2) and c(q1) = c(q2) and d(q1) = d(q2);

    % writes a quaternion                                                    %
    procedure writeonQ( reference(Quaternion) value q ) ;
        writeon( "(", a(q), ", ", b(q), ", ", c(q), ", ", d(q), ")" );


    % test q1q2 = q2q1                                                       %
    reference(Quaternion) q, q1, q2;

    q  := Quaternion( 1, 2, 3, 4 );
    q1 := Quaternion( 2, 3, 4, 5 );
    q2 := Quaternion( 3, 4, 5, 6 );

    % set output format                                                      %
    s_w := 0; r_format := "A"; r_w := 5; r_d := 1;

    write( "      q:" );writeonQ( q );
    write( "     q1:" );writeonQ( q1 );
    write( "     q2:" );writeonQ( q2 );
    write( "norm  q:" );writeon( normQ( q ) );
    write( "norm q1:" );writeon( normQ( q1 ) );
    write( "norm q2:" );writeon( normQ( q2 ) );

    write( " conj q:" );writeonQ( conjQ( q ) );
    write( "    - q:" );writeonQ( negQ( q ) );
    write( "  7 + q:" );writeonQ( addRQ( 7, q ) );
    write( "  q + 9:" );writeonQ( addQR( q, 9 ) );
    write( " q + q1:" );writeonQ( addQQ( q, q1 ) );
    write( "  3 * q:" );writeonQ( mulRQ( 3, q ) );
    write( "  q * 4:" );writeonQ( mulQR( q, 4 ) );

    % check that q1q2 not = q2q1                                             %
    if equalQ( mulQQ( q1, q2 ), mulQQ( q2, q1 ) )
    then write( "q1q2 = q2q1 ??" )
    else write( "q1q2 <> q2q1" );

    write( "   q1q2:" );writeonQ( mulQQ( q1, q2 ) );
    write( "   q2q1:" );writeonQ( mulQQ( q2, q1 ) );

end.
Output:
      q:(  1.0,   2.0,   3.0,   4.0)
     q1:(  2.0,   3.0,   4.0,   5.0)
     q2:(  3.0,   4.0,   5.0,   6.0)
norm  q:  5.4
norm q1:  7.3
norm q2:  9.2
 conj q:(  1.0,  -2.0,  -3.0,  -4.0)
    - q:( -1.0,  -2.0,  -3.0,  -4.0)
  7 + q:(  8.0,   2.0,   3.0,   4.0)
  q + 9:( 10.0,   2.0,   3.0,   4.0)
 q + q1:(  3.0,   5.0,   7.0,   9.0)
  3 * q:(  3.0,   6.0,   9.0,  12.0)
  q * 4:(  4.0,   8.0,  12.0,  16.0)
q1q2 <> q2q1
   q1q2:(-56.0,  16.0,  24.0,  26.0)
   q2q1:(-56.0,  18.0,  20.0,  28.0)

Arturo

qnorm: $ => [sqrt fold & [x y] -> x + y*y]

qneg: $ => [map & => neg]

qconj: $[q] [@[q\0] ++ qneg drop q]

qaddr: function [q r][
    [a b c d]: q
    @[a+r b c d]
]

qadd: $ => [map couple & & => sum]

qmulr: $[q r] [map q'x -> x*r]

qmul: function [q1 q2][
    [a1 b1 c1 d1]: q1
    [a2 b2 c2 d2]: q2
    @[
        (((a1*a2) - b1*b2) - c1*c2) - d1*d2,
        (((a1*b2) + b1*a2) + c1*d2) - d1*c2,
        (((a1*c2) - b1*d2) + c1*a2) + d1*b2,
        (((a1*d2) + b1*c2) - c1*b2) + d1*a2
    ]
]

; --- test quaternions ---
q:  [1 2 3 4]
q1: [2 3 4 5]
q2: [3 4 5 6]
r: 7

print ['qnorm q '= qnorm q]
print ['qneg q '= qneg q]
print ['qconj q '= qconj q]
print ['qaddr q r '= qaddr q r]
print ['qmulr q r '= qmulr q r]
print ['qadd q1 q2 '= qadd q1 q2]
print ['qmul q1 q2 '= qmul q1 q2]
print ['qmul q2 q1 '= qmul q2 q1]
Output:
qnorm [1 2 3 4] = 5.477225575051661 
qneg [1 2 3 4] = [-1 -2 -3 -4] 
qconj [1 2 3 4] = [1 -2 -3 -4] 
qaddr [1 2 3 4] 7 = [8 2 3 4] 
qmulr [1 2 3 4] 7 = [7 14 21 28] 
qadd [2 3 4 5] [3 4 5 6] = [5 7 9 11] 
qmul [2 3 4 5] [3 4 5 6] = [-56 16 24 26] 
qmul [3 4 5 6] [2 3 4 5] = [-56 18 20 28]

ATS

Library: ats2-xprelude
//--------------------------------------------------------------------

#include "share/atspre_staload.hats"

//--------------------------------------------------------------------

(* Here is one way to get a sqrt function without going beyond the ATS
   prelude. The prelude (at the time of this writing) contains some
   templates for which implementations were never added. Here I add an
   implementation.

   The ats2-xprelude package at
   https://sourceforge.net/p/chemoelectric/ats2-xprelude contains a
   much more extensive and natural interface to the C math library. *)

%{^
#include <math.h>
%}

implement                       (* "Generic" square root. *)
gsqrt_val<double> x =
  (* Call "sqrt" from the C math library. *)
  $extfcall (double, "sqrt", x)

//--------------------------------------------------------------------

abst@ype quaternion (tk : tkind) =
  (* The following determines the SIZE of a quaternion, but not its
     actual representation: *)
  @(g0float tk, g0float tk, g0float tk, g0float tk)

extern fn {tk : tkind} quaternion_make :
  (g0float tk, g0float tk, g0float tk, g0float tk) -<> quaternion tk

extern fn {tk : tkind} fprint_quaternion :
  (FILEref, quaternion tk) -> void
extern fn {tk : tkind} print_quaternion :
  quaternion tk -> void

extern fn {tk : tkind} quaternion_norm_squared :
  quaternion tk -<> g0float tk
extern fn {tk : tkind} quaternion_norm :
  quaternion tk -< !exn > g0float tk

extern fn {tk : tkind} quaternion_neg :
  quaternion tk -<> quaternion tk
extern fn {tk : tkind} quaternion_conj :
  quaternion tk -<> quaternion tk

extern fn {tk : tkind} add_quaternion_g0float :
  (quaternion tk, g0float tk) -<> quaternion tk
extern fn {tk : tkind} add_g0float_quaternion :
  (g0float tk, quaternion tk) -<> quaternion tk
extern fn {tk : tkind} add_quaternion_quaternion :
  (quaternion tk, quaternion tk) -<> quaternion tk

extern fn {tk : tkind} mul_quaternion_g0float :
  (quaternion tk, g0float tk) -<> quaternion tk
extern fn {tk : tkind} mul_g0float_quaternion :
  (g0float tk, quaternion tk) -<> quaternion tk
extern fn {tk : tkind} mul_quaternion_quaternion :
  (quaternion tk, quaternion tk) -<> quaternion tk

extern fn {tk : tkind} quaternion_eq :
  (quaternion tk, quaternion tk) -<> bool

overload fprint with fprint_quaternion
overload print with print_quaternion

overload norm_squared with quaternion_norm_squared
overload norm with quaternion_norm

overload ~ with quaternion_neg
overload conj with quaternion_conj

overload + with add_quaternion_g0float
overload + with add_g0float_quaternion
overload + with add_quaternion_quaternion

overload * with mul_quaternion_g0float
overload * with mul_g0float_quaternion
overload * with mul_quaternion_quaternion

overload = with quaternion_eq

//--------------------------------------------------------------------

local

  (* Now we decide the REPRESENTATION of a quaternion. A quaternion is
     represented as an unboxed 4-tuple of "real" numbers of any one
     particular typekind. *)
  typedef _quaternion (tk : tkind) =
    @(g0float tk, g0float tk, g0float tk, g0float tk)

  assume quaternion tk = _quaternion tk

in (* local *)

  implement {tk}
  quaternion_make (a, b, c, d) =
    @(a, b, c, d)  

  implement {tk}
  fprint_quaternion (outf, q) =
    let
      typedef t = g0float tk
      val @(a, b, c, d) = q
    in
      fprint_val<t> (outf, a);
      if g0i2f 0 <= b then fprint_val<string> (outf, "+");
      fprint_val<t> (outf, b);
      fprint_val<string> (outf, "i");
      if g0i2f 0 <= c then fprint_val<string> (outf, "+");
      fprint_val<t> (outf, c);
      fprint_val<string> (outf, "j");
      if g0i2f 0 <= d then fprint_val<string> (outf, "+");
      fprint_val<t> (outf, d);
      fprint_val<string> (outf, "k");
    end

  implement {tk}
  print_quaternion q =
    fprint_quaternion (stdout_ref, q)

  implement {tk}
  quaternion_norm_squared q =
    let
      val @(a, b, c, d) = q
    in
      (a * a) + (b * b) + (c * c) + (d * d)
    end

  implement {tk}
  quaternion_norm q =
    gsqrt_val<g0float tk> (quaternion_norm_squared q)

  implement {tk}
  quaternion_neg q =
    let
      val @(a, b, c, d) = q
    in
      @(~a, ~b, ~c, ~d)
    end    

  implement {tk}
  quaternion_conj q =
    let
      val @(a, b, c, d) = q
    in
      @(a, ~b, ~c, ~d)
    end    

  implement {tk}
  add_quaternion_g0float (q, r) =
    let
      val @(a, b, c, d) = q
    in
      @(a + r, b, c, d)
    end    

  implement {tk}
  add_g0float_quaternion (r, q) =
    let
      val @(a, b, c, d) = q
    in
      @(r + a, b, c, d)
    end    

  implement {tk}
  add_quaternion_quaternion (q1, q2) =
    let
      val @(a1, b1, c1, d1) = q1
      and @(a2, b2, c2, d2) = q2
    in
      @(a1 + a2, b1 + b2, c1 + c2, d1 + d2)
    end    

  implement {tk}
  mul_quaternion_g0float (q, r) =
    let
      val @(a, b, c, d) = q
    in
      @(a * r, b * r, c * r, d * r)
    end    

  implement {tk}
  mul_g0float_quaternion (r, q) =
    let
      val @(a, b, c, d) = q
    in
      @(r * a, r * b, r * c, r * d)
    end    

  implement {tk}
  mul_quaternion_quaternion (q1, q2) =
    let
      val @(a1, b1, c1, d1) = q1
      and @(a2, b2, c2, d2) = q2
    in
      @((a1 * a2) - (b1 * b2) - (c1 * c2) - (d1 * d2),
        (a1 * b2) + (b1 * a2) + (c1 * d2) - (d1 * c2),
        (a1 * c2) - (b1 * d2) + (c1 * a2) + (d1 * b2),
        (a1 * d2) + (b1 * c2) - (c1 * b2) + (d1 * a2))
    end    

  implement {tk}
  quaternion_eq (q1, q2) =
    let
      val @(a1, b1, c1, d1) = q1
      and @(a2, b2, c2, d2) = q2
    in
      (a1 = a2) * (b1 = b2) * (c1 = c2) * (d1 = d2)
    end

end (* local *)

//--------------------------------------------------------------------

val q = quaternion_make (1.0, 2.0, 3.0, 4.0)
and q1 = quaternion_make (2.0, 3.0, 4.0, 5.0)
and q2 = quaternion_make (3.0, 4.0, 5.0, 6.0)
and r = 7.0

implement
main0 () =
  let
    (* Let us print double precision numbers in a format more readable
       than is the prelude's default. *)
    implement
    fprint_val<double> (outf, x) =
      let
        typedef f = $extype"FILE *"
        val _ = $extfcall (int, "fprintf", $UNSAFE.cast{f} outf,
                           "%g", x)
      in
      end
  in
    println! ("q = ", q);
    println! ("q1 = ", q1);
    println! ("q2 = ", q2);
    println! ();
    println! ("||q|| = ", norm q);
    println! ("||q1|| = ", norm q1);
    println! ("||q2|| = ", norm q2);
    println! ();
    println! ("-q = ", ~q);
    println! ("-q1 = ", ~q1);
    println! ("-q2 = ", ~q2);
    println! ();
    println! ("conj q = ", conj q);
    println! ("conj q1 = ", conj q1);
    println! ("conj q2 = ", conj q2);
    println! ();
    println! ("q + r = ", q + r);
    println! ("r + q = ", r + q);
    println! ("q1 + q2 = ", q1 + q2);
    println! ();
    println! ("q * r = ", q * r);
    println! ("r * q = ", r * q);
    println! ("q1 * q2 = ", q1 * q2);
    println! ("q2 * q1 = ", q2 * q1);
    println! ("((q1 * q2) = (q2 * q1)) is ", (q1 * q2) = (q2 * q1))
  end

//--------------------------------------------------------------------
Output:
$ patscc -std=gnu2x -O2 quaternions_task.dats -lm && ./a.out
q = 1+2i+3j+4k
q1 = 2+3i+4j+5k
q2 = 3+4i+5j+6k

||q|| = 5.477226
||q1|| = 7.348469
||q2|| = 9.273618

-q = -1-2i-3j-4k
-q1 = -2-3i-4j-5k
-q2 = -3-4i-5j-6k

conj q = 1-2i-3j-4k
conj q1 = 2-3i-4j-5k
conj q2 = 3-4i-5j-6k

q + r = 8+2i+3j+4k
r + q = 8+2i+3j+4k
q1 + q2 = 5+7i+9j+11k

q * r = 7+14i+21j+28k
r * q = 7+14i+21j+28k
q1 * q2 = -56+16i+24j+26k
q2 * q1 = -56+18i+20j+28k
((q1 * q2) = (q2 * q1)) is false

AutoHotkey

Works with: AutoHotkey_L

(AutoHotkey1.1+)

q  := [1, 2, 3, 4]
q1 := [2, 3, 4, 5]
q2 := [3, 4, 5, 6]
r := 7

MsgBox, % "q = " PrintQ(q)
	. "`nq1 = " PrintQ(q1)
	. "`nq2 = " PrintQ(q2)
	. "`nr = " r
	. "`nNorm(q) = " Norm(q)
	. "`nNegative(q) = " PrintQ(Negative(q))
	. "`nConjugate(q) = " PrintQ(Conjugate(q))
	. "`nq + r = " PrintQ(AddR(q, r))
	. "`nq1 + q2 = " PrintQ(AddQ(q1, q2))
	. "`nq2 + q1 = " PrintQ(AddQ(q2, q1))
	. "`nqr = " PrintQ(MulR(q, r))
	. "`nq1q2 = " PrintQ(MulQ(q1, q2))
	. "`nq2q1 = " PrintQ(MulQ(q2, q1))

Norm(q) {
	return sqrt(q[1]**2 + q[2]**2 + q[3]**2 + q[4]**2)
}

Negative(q) {
	a := []
	for k, v in q
		a[A_Index] := v * -1
	return a
}

Conjugate(q) {
	a := []
	for k, v in q
		a[A_Index] := v * (A_Index = 1 ? 1 : -1)
	return a
}

AddR(q, r) {
	a := []
	for k, v in q 
		a[A_Index] := v + (A_Index = 1 ? r : 0)
	return a
}

AddQ(q1, q2) {
	a := []
	for k, v in q1
		a[A_Index] := v + q2[A_Index]
	return a
}

MulR(q, r) {
	a := []
	for k, v in q
		a[A_Index] := v * r
	return a
}

MulQ(q, u) {
	a := []
	, a[1] := q[1]*u[1] - q[2]*u[2] - q[3]*u[3] - q[4]*u[4]
	, a[2] := q[1]*u[2] + q[2]*u[1] + q[3]*u[4] - q[4]*u[3]
	, a[3] := q[1]*u[3] - q[2]*u[4] + q[3]*u[1] + q[4]*u[2]
	, a[4] := q[1]*u[4] + q[2]*u[3] - q[3]*u[2] + q[4]*u[1]
	return a
}

PrintQ(q, b="(") {
	for k, v in q
		b .= v (A_Index = q.MaxIndex() ? ")" : ", ")
	return b
}
Output:
q = (1, 2, 3, 4)
q1 = (2, 3, 4, 5)
q2 = (3, 4, 5, 6)
r = 7
Norm(q) = 5.477226
Negative(q) = (-1, -2, -3, -4)
Conjugate(q) = (1, -2, -3, -4)
q + r = (8, 2, 3, 4)
q1 + q2 = (5, 7, 9, 11)
q2 + q1 = (5, 7, 9, 11)
qr = (7, 14, 21, 28)
q1q2 = (-56, 16, 24, 26)
q2q1 = (-56, 18, 20, 28)

Axiom

Axiom has built-in support for quaternions.

qi := quatern$Quaternion(Integer);

             Type: ((Integer,Integer,Integer,Integer) -> Quaternion(Integer))
q  := qi(1,2,3,4);

                                                    Type: Quaternion(Integer)
q1 := qi(2,3,4,5);

                                                    Type: Quaternion(Integer)
q2 := qi(3,4,5,6);

                                                    Type: Quaternion(Integer)
r : Integer := 7;

                                                                Type: Integer
sqrt norm q

         +--+
   (6)  \|30
                                                        Type: AlgebraicNumber
-q

   (7)  - 1 - 2i - 3j - 4k
                                                    Type: Quaternion(Integer)
conjugate q

   (8)  1 - 2i - 3j - 4k
                                                    Type: Quaternion(Integer)
r + q

   (9)  8 + 2i + 3j + 4k
                                                    Type: Quaternion(Integer)
q1 + q2

   (10)  5 + 7i + 9j + 11k
                                                    Type: Quaternion(Integer)
q*r

   (11)  7 + 14i + 21j + 28k
                                                    Type: Quaternion(Integer)
r*q

   (12)  7 + 14i + 21j + 28k
                                                    Type: Quaternion(Integer)
q1*q2 ~= q2*q1

   (13)  true
                                                                Type: Boolean

BASIC

BASIC256

Works with: BASIC256 version 2.0.0.11
dim q(4)
dim q1(4)
dim q2(4)
q[0] = 1: q[1] = 2: q[2] = 3: q[3] = 4
q1[0] = 2: q1[1] = 3: q1[2] = 4: q1[3] = 5
q2[0] = 3: q2[1] = 4: q2[2] = 5: q2[3] = 6
r = 7

function printq(q)
 return "("+q[0]+", "+q[1]+", "+q[2]+", "+q[3]+")"
end function

function q_equal(q1, q2)
 return q1[0]=q2[0] and q1[1]=q2[1] and q1[2]=q2[2] and q1[3]=q2[3]
end function

function q_norm(q)
 return sqr(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3])
end function

function q_neg(q)
 dim result[4]
 result[0] = -q[0]
 result[1] = -q[1]
 result[2] = -q[2]
 result[3] = -q[3]
 return result
end function

function q_conj(q)
 dim result[4]
 result[0] = q[0]
 result[1] = -q[1]
 result[2] = -q[2]
 result[3] = -q[3]
 return result
end function

function q_addreal(q, r)
 dim result[4]
 result[0] = q[0]+r
 result[1] = q[1]
 result[2] = q[2]
 result[3] = q[3]
 return result
end function

function q_add(q1, q2)
 dim result[4]
 result[0] = q1[0]+q2[0]
 result[1] = q1[1]+q2[1]
 result[2] = q1[2]+q2[2]
 result[3] = q1[3]+q2[3]
 return result
end function

function q_mulreal(q, r)
 dim result[4]
 result[0] = q[0]*r
 result[1] = q[1]*r
 result[2] = q[2]*r
 result[3] = q[3]*r
 return result
end function

function q_mul(q1, q2)
 dim result[4]
 result[0] = q1[0]*q2[0]-q1[1]*q2[1]-q1[2]*q2[2]-q1[3]*q2[3]
 result[1] = q1[0]*q2[1]+q1[1]*q2[0]+q1[2]*q2[3]-q1[3]*q2[2]
 result[2] = q1[0]*q2[2]-q1[1]*q2[3]+q1[2]*q2[0]+q1[3]*q2[1]
 result[3] = q1[0]*q2[3]+q1[1]*q2[2]-q1[2]*q2[1]+q1[3]*q2[0]
 return result
end function

print "q = ";printq(q)
print "q1 = ";printq(q1)
print "q2 = ";printq(q2)
print "r = "; r
print "norm(q) = "; q_norm(q)
print "neg(q) = ";printq(q_neg(q))
print "conjugate(q) = ";printq(q_conj(q))
print "q+r = ";printq(q_addreal(q,r))
print "q1+q2 = ";printq(q_add(q1,q2))
print "qr = ";printq(q_mulreal(q,r))
print "q1q2 = ";printq(q_mul(q1,q2))
print "q2q1 = ";printq(q_mul(q2,q1))
Output:
q = (1, 2, 3, 4)
q1 = (2, 3, 4, 5)
q2 = (3, 4, 5, 6)
r = 7
norm(q) = 5.47722557505
neg(q) = (-1, -2, -3, -4)
conjugate(q) = (1, -2, -3, -4)
q+r = (8, 2, 3, 4)
q1+q2 = (5, 7, 9, 11)
qr = (7, 14, 21, 28)
q1q2 = (-56, 16, 24, 26)
q2q1 = (-56, 18, 20, 28)

BBC BASIC

Although BBC BASIC doesn't have native support for quaternions its array arithmetic provides all of the required operations either directly or very straightforwardly.

      DIM q(3), q1(3), q2(3), t(3)
      q() = 1, 2, 3, 4
      q1() = 2, 3, 4, 5
      q2() = 3, 4, 5, 6
      r = 7
      
      PRINT "q = " FNq_show(q())
      PRINT "q1 = " FNq_show(q1())
      PRINT "q2 = " FNq_show(q2())
      PRINT "r = "; r
      PRINT "norm(q) = "; FNq_norm(q())
      t() = q() : PROCq_neg(t()) : PRINT "neg(q) = " FNq_show(t())
      t() = q() : PROCq_conj(t()) : PRINT "conjugate(q) = " FNq_show(t())
      t() = q() : PROCq_addreal(t(),r) : PRINT "q + r = " FNq_show(t())
      t() = q1() : PROCq_add(t(),q2()) : PRINT "q1 + q2 = " FNq_show(t())
      t() = q2() : PROCq_add(t(),q1()) : PRINT "q2 + q1 = " FNq_show(t())
      t() = q() : PROCq_mulreal(t(),r) : PRINT "qr = " FNq_show(t())
      t() = q1() : PROCq_mul(t(),q2()) : PRINT "q1q2 = " FNq_show(t())
      t() = q2() : PROCq_mul(t(),q1()) : PRINT "q2q1 = " FNq_show(t())
      END
      
      DEF FNq_norm(q()) = MOD(q())
      
      DEF PROCq_neg(q()) : q() *= -1 : ENDPROC
      
      DEF PROCq_conj(q()) : q() *= -1 : q(0) *= -1 : ENDPROC
      
      DEF PROCq_addreal(q(), r) : q(0) += r : ENDPROC
      
      DEF PROCq_add(q(), r()) : q() += r() : ENDPROC
      
      DEF PROCq_mulreal(q(), r) : q() *= r : ENDPROC
      
      DEF PROCq_mul(q(), r()) : LOCAL s() : DIM s(3,3)
      s() = r(0), -r(1), -r(2), -r(3), r(1), r(0),  r(3), -r(2), \
      \     r(2), -r(3),  r(0),  r(1), r(3), r(2), -r(1),  r(0)
      q() = s() . q()
      ENDPROC
      
      DEF FNq_show(q()) : LOCAL i%, a$ : a$ = "("
      FOR i% = 0 TO 3 : a$ += STR$(q(i%)) + ", " : NEXT
      = LEFT$(LEFT$(a$)) + ")"
Output:
q = (1, 2, 3, 4)
q1 = (2, 3, 4, 5)
q2 = (3, 4, 5, 6)
r = 7
norm(q) = 5.47722558
neg(q) = (-1, -2, -3, -4)
conjugate(q) = (1, -2, -3, -4)
q + r = (8, 2, 3, 4)
q1 + q2 = (5, 7, 9, 11)
q2 + q1 = (5, 7, 9, 11)
qr = (7, 14, 21, 28)
q1q2 = (-56, 16, 24, 26)
q2q1 = (-56, 18, 20, 28)

C

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>

typedef struct quaternion
{
  double q[4];
} quaternion_t;


quaternion_t *quaternion_new(void)
{
  return malloc(sizeof(quaternion_t));
}

quaternion_t *quaternion_new_set(double q1,
				 double q2,
				 double q3,
				 double q4)
{
  quaternion_t *q = malloc(sizeof(quaternion_t));
  if (q != NULL) {
    q->q[0] = q1; q->q[1] = q2; q->q[2] = q3; q->q[3] = q4;
  }
  return q;
}


void quaternion_copy(quaternion_t *r, quaternion_t *q)
{
  size_t i;

  if (r == NULL || q == NULL) return;
  for(i = 0; i < 4; i++) r->q[i] = q->q[i];
}


double quaternion_norm(quaternion_t *q)
{
  size_t i;
  double r = 0.0;
  
  if (q == NULL) {
    fprintf(stderr, "NULL quaternion in norm\n");
    return 0.0;
  }

  for(i = 0; i < 4; i++) r += q->q[i] * q->q[i];
  return sqrt(r);
}


void quaternion_neg(quaternion_t *r, quaternion_t *q)
{
  size_t i;

  if (q == NULL || r == NULL) return;
  for(i = 0; i < 4; i++) r->q[i] = -q->q[i];
}


void quaternion_conj(quaternion_t *r, quaternion_t *q)
{
  size_t i;

  if (q == NULL || r == NULL) return;
  r->q[0] = q->q[0];
  for(i = 1; i < 4; i++) r->q[i] = -q->q[i];
}


void quaternion_add_d(quaternion_t *r, quaternion_t *q, double d)
{
  if (q == NULL || r == NULL) return;
  quaternion_copy(r, q);
  r->q[0] += d;
}


void quaternion_add(quaternion_t *r, quaternion_t *a, quaternion_t *b)
{
  size_t i;

  if (r == NULL || a == NULL || b == NULL) return;
  for(i = 0; i < 4; i++) r->q[i] = a->q[i] + b->q[i];
}


void quaternion_mul_d(quaternion_t *r, quaternion_t *q, double d)
{
  size_t i;

  if (r == NULL || q == NULL) return;
  for(i = 0; i < 4; i++) r->q[i] = q->q[i] * d;
}

bool quaternion_equal(quaternion_t *a, quaternion_t *b)
{
  size_t i;
  
  for(i = 0; i < 4; i++) if (a->q[i] != b->q[i]) return false;
  return true;
}


#define A(N) (a->q[(N)])
#define B(N) (b->q[(N)])
#define R(N) (r->q[(N)])
void quaternion_mul(quaternion_t *r, quaternion_t *a, quaternion_t *b)
{
  size_t i;
  double ri = 0.0;

  if (r == NULL || a == NULL || b == NULL) return;
  R(0) = A(0)*B(0) - A(1)*B(1) - A(2)*B(2) - A(3)*B(3);
  R(1) = A(0)*B(1) + A(1)*B(0) + A(2)*B(3) - A(3)*B(2);
  R(2) = A(0)*B(2) - A(1)*B(3) + A(2)*B(0) + A(3)*B(1);
  R(3) = A(0)*B(3) + A(1)*B(2) - A(2)*B(1) + A(3)*B(0);
}
#undef A
#undef B
#undef R


void quaternion_print(quaternion_t *q)
{
  if (q == NULL) return;
  printf("(%lf, %lf, %lf, %lf)\n", 
	 q->q[0], q->q[1], q->q[2], q->q[3]);
}
int main()
{
  size_t i;
  double d = 7.0;
  quaternion_t *q[3];
  quaternion_t *r  = quaternion_new();

  quaternion_t *qd = quaternion_new_set(7.0, 0.0, 0.0, 0.0);
  q[0] = quaternion_new_set(1.0, 2.0, 3.0, 4.0);
  q[1] = quaternion_new_set(2.0, 3.0, 4.0, 5.0);
  q[2] = quaternion_new_set(3.0, 4.0, 5.0, 6.0);

  printf("r = %lf\n", d);
  
  for(i = 0; i < 3; i++) {
    printf("q[%u] = ", i);
    quaternion_print(q[i]);
    printf("abs q[%u] = %lf\n", i, quaternion_norm(q[i]));
  }

  printf("-q[0] = ");
  quaternion_neg(r, q[0]);
  quaternion_print(r);

  printf("conj q[0] = ");
  quaternion_conj(r, q[0]);
  quaternion_print(r);

  printf("q[1] + q[2] = ");
  quaternion_add(r, q[1], q[2]);
  quaternion_print(r);

  printf("q[2] + q[1] = ");
  quaternion_add(r, q[2], q[1]);
  quaternion_print(r);
  

  printf("q[0] * r = ");
  quaternion_mul_d(r, q[0], d);
  quaternion_print(r);

  printf("q[0] * (r, 0, 0, 0) = ");
  quaternion_mul(r, q[0], qd);
  quaternion_print(r);
  

  printf("q[1] * q[2] = ");
  quaternion_mul(r, q[1], q[2]);
  quaternion_print(r);

  printf("q[2] * q[1] = ");
  quaternion_mul(r, q[2], q[1]);
  quaternion_print(r);  


  free(q[0]); free(q[1]); free(q[2]); free(r);
  return EXIT_SUCCESS;
}

C#

using System;

struct Quaternion : IEquatable<Quaternion>
{
    public readonly double A, B, C, D;

    public Quaternion(double a, double b, double c, double d)
    {
        this.A = a;
        this.B = b;
        this.C = c;
        this.D = d;
    }

    public double Norm()
    {
        return Math.Sqrt(A * A + B * B + C * C + D * D);
    }

    public static Quaternion operator -(Quaternion q)
    {
        return new Quaternion(-q.A, -q.B, -q.C, -q.D);
    }

    public Quaternion Conjugate()
    {
        return new Quaternion(A, -B, -C, -D);
    }

    // implicit conversion takes care of real*quaternion and real+quaternion
    public static implicit operator Quaternion(double d)
    {
        return new Quaternion(d, 0, 0, 0);
    }

    public static Quaternion operator +(Quaternion q1, Quaternion q2)
    {
        return new Quaternion(q1.A + q2.A, q1.B + q2.B, q1.C + q2.C, q1.D + q2.D);
    }

    public static Quaternion operator *(Quaternion q1, Quaternion q2)
    {
        return new Quaternion(
            q1.A * q2.A - q1.B * q2.B - q1.C * q2.C - q1.D * q2.D,
            q1.A * q2.B + q1.B * q2.A + q1.C * q2.D - q1.D * q2.C,
            q1.A * q2.C - q1.B * q2.D + q1.C * q2.A + q1.D * q2.B,
            q1.A * q2.D + q1.B * q2.C - q1.C * q2.B + q1.D * q2.A);
    }

    public static bool operator ==(Quaternion q1, Quaternion q2)
    {
        return q1.A == q2.A && q1.B == q2.B && q1.C == q2.C && q1.D == q2.D;
    }

    public static bool operator !=(Quaternion q1, Quaternion q2)
    {
        return !(q1 == q2);
    }

    #region Object Members

    public override bool Equals(object obj)
    {
        if (obj is Quaternion)
            return Equals((Quaternion)obj);

        return false;
    }

    public override int GetHashCode()
    {
        return A.GetHashCode() ^ B.GetHashCode() ^ C.GetHashCode() ^ D.GetHashCode();
    }

    public override string ToString()
    {
        return string.Format("Q({0}, {1}, {2}, {3})", A, B, C, D);
    }

    #endregion

    #region IEquatable<Quaternion> Members

    public bool Equals(Quaternion other)
    {
        return other == this;
    }

    #endregion
}

Demonstration:

using System;

static class Program
{
    static void Main(string[] args)
    {
        Quaternion q = new Quaternion(1, 2, 3, 4);
        Quaternion q1 = new Quaternion(2, 3, 4, 5);
        Quaternion q2 = new Quaternion(3, 4, 5, 6);
        double r = 7;

        Console.WriteLine("q = {0}", q);
        Console.WriteLine("q1 = {0}", q1);
        Console.WriteLine("q2 = {0}", q2);
        Console.WriteLine("r = {0}", r);

        Console.WriteLine("q.Norm() = {0}", q.Norm());
        Console.WriteLine("q1.Norm() = {0}", q1.Norm());
        Console.WriteLine("q2.Norm() = {0}", q2.Norm());

        Console.WriteLine("-q = {0}", -q);
        Console.WriteLine("q.Conjugate() = {0}", q.Conjugate());

        Console.WriteLine("q + r = {0}", q + r);
        Console.WriteLine("q1 + q2 = {0}", q1 + q2);
        Console.WriteLine("q2 + q1 = {0}", q2 + q1);

        Console.WriteLine("q * r = {0}", q * r);
        Console.WriteLine("q1 * q2 = {0}", q1 * q2);
        Console.WriteLine("q2 * q1 = {0}", q2 * q1);

        Console.WriteLine("q1*q2 {0} q2*q1", (q1 * q2) == (q2 * q1) ? "==" : "!=");
    }
}
Output:
q = Q(1, 2, 3, 4)
q1 = Q(2, 3, 4, 5)
q2 = Q(3, 4, 5, 6)
r = 7
q.Norm() = 5.47722557505166
q1.Norm() = 7.34846922834953
q2.Norm() = 9.2736184954957
-q = Q(-1, -2, -3, -4)
q.Conjugate() = Q(1, -2, -3, -4)
q + r = Q(8, 2, 3, 4)
q1 + q2 = Q(5, 7, 9, 11)
q2 + q1 = Q(5, 7, 9, 11)
q * r = Q(7, 14, 21, 28)
q1 * q2 = Q(-56, 16, 24, 26)
q2 * q1 = Q(-56, 18, 20, 28)
q1*q2 != q2*q1

C++

This example uses templates to provide the underlying data-type, and includes several extra functions and constructors that often come up when using quaternions.

#include <iostream>
using namespace std;

template<class T = double>
class Quaternion
{
public:
  T w, x, y, z;

  // Numerical constructor
  Quaternion(const T &w, const T &x, const T &y, const T &z): w(w), x(x), y(y), z(z) {};
  Quaternion(const T &x, const T &y, const T &z): w(T()), x(x), y(y), z(z) {}; // For 3-rotations
  Quaternion(const T &r): w(r), x(T()), y(T()), z(T()) {};
  Quaternion(): w(T()), x(T()), y(T()), z(T()) {};

  // Copy constructor and assignment
  Quaternion(const Quaternion &q): w(q.w), x(q.x), y(q.y), z(q.z) {};
  Quaternion& operator=(const Quaternion &q) { w=q.w; x=q.x; y=q.y; z=q.z; return *this; }

  // Unary operators
  Quaternion operator-() const { return Quaternion(-w, -x, -y, -z); }
  Quaternion operator~() const { return Quaternion(w, -x, -y, -z); } // Conjugate

  // Norm-squared. SQRT would have to be made generic to be used here
  T normSquared() const { return w*w + x*x + y*y + z*z; }

  // In-place operators
  Quaternion& operator+=(const T &r) 
    { w += r; return *this; }
  Quaternion& operator+=(const Quaternion &q) 
    { w += q.w; x += q.x; y += q.y; z += q.z; return *this; }

  Quaternion& operator-=(const T &r) 
    { w -= r; return *this; }
  Quaternion& operator-=(const Quaternion &q) 
    { w -= q.w; x -= q.x; y -= q.y; z -= q.z; return *this; }

  Quaternion& operator*=(const T &r) 
    { w *= r; x *= r; y *= r; z *= r; return *this; }
  Quaternion& operator*=(const Quaternion &q) 
  { 
    T oldW(w), oldX(x), oldY(y), oldZ(z);
    w = oldW*q.w - oldX*q.x - oldY*q.y - oldZ*q.z;
    x = oldW*q.x + oldX*q.w + oldY*q.z - oldZ*q.y;
    y = oldW*q.y + oldY*q.w + oldZ*q.x - oldX*q.z;
    z = oldW*q.z + oldZ*q.w + oldX*q.y - oldY*q.x;
    return *this;
  }
    
  Quaternion& operator/=(const T &r) 
    { w /= r; x /= r; y /= r; z /= r; return *this; }
  Quaternion& operator/=(const Quaternion &q) 
  { 
    T oldW(w), oldX(x), oldY(y), oldZ(z), n(q.normSquared());
    w = (oldW*q.w + oldX*q.x + oldY*q.y + oldZ*q.z) / n;
    x = (oldX*q.w - oldW*q.x + oldY*q.z - oldZ*q.y) / n;
    y = (oldY*q.w - oldW*q.y + oldZ*q.x - oldX*q.z) / n;
    z = (oldZ*q.w - oldW*q.z + oldX*q.y - oldY*q.x) / n;
    return *this;
  }

  // Binary operators based on in-place operators
  Quaternion operator+(const T &r) const { return Quaternion(*this) += r; }
  Quaternion operator+(const Quaternion &q) const { return Quaternion(*this) += q; }
  Quaternion operator-(const T &r) const { return Quaternion(*this) -= r; }
  Quaternion operator-(const Quaternion &q) const { return Quaternion(*this) -= q; }
  Quaternion operator*(const T &r) const { return Quaternion(*this) *= r; }
  Quaternion operator*(const Quaternion &q) const { return Quaternion(*this) *= q; }
  Quaternion operator/(const T &r) const { return Quaternion(*this) /= r; }
  Quaternion operator/(const Quaternion &q) const { return Quaternion(*this) /= q; }

  // Comparison operators, as much as they make sense
  bool operator==(const Quaternion &q) const 
    { return (w == q.w) && (x == q.x) && (y == q.y) && (z == q.z); }
  bool operator!=(const Quaternion &q) const { return !operator==(q); }

  // The operators above allow quaternion op real. These allow real op quaternion.
  // Uses the above where appropriate.
  template<class T> friend Quaternion<T> operator+(const T &r, const Quaternion<T> &q);
  template<class T> friend Quaternion<T> operator-(const T &r, const Quaternion<T> &q);
  template<class T> friend Quaternion<T> operator*(const T &r, const Quaternion<T> &q);
  template<class T> friend Quaternion<T> operator/(const T &r, const Quaternion<T> &q);
  
  // Allows cout << q 
  template<class T> friend ostream& operator<<(ostream &io, const Quaternion<T> &q);
};

// Friend functions need to be outside the actual class definition
template<class T>
Quaternion<T> operator+(const T &r, const Quaternion<T> &q) 
  { return q+r; }

template<class T>
Quaternion<T> operator-(const T &r, const Quaternion<T> &q)
  { return Quaternion<T>(r-q.w, q.x, q.y, q.z); }

template<class T>
Quaternion<T> operator*(const T &r, const Quaternion<T> &q) 
  { return q*r; }

template<class T>
Quaternion<T> operator/(const T &r, const Quaternion<T> &q)
{
  T n(q.normSquared());
  return Quaternion(r*q.w/n, -r*q.x/n, -r*q.y/n, -r*q.z/n);
}

template<class T>
ostream& operator<<(ostream &io, const Quaternion<T> &q)
{ 
  io << q.w;
  (q.x < T()) ? (io << " - " << (-q.x) << "i") : (io << " + " << q.x << "i");
  (q.y < T()) ? (io << " - " << (-q.y) << "j") : (io << " + " << q.y << "j");
  (q.z < T()) ? (io << " - " << (-q.z) << "k") : (io << " + " << q.z << "k");
  return io;
}

Test program:

int main()
{
  Quaternion<> q0(1, 2, 3, 4);
  Quaternion<> q1(2, 3, 4, 5);
  Quaternion<> q2(3, 4, 5, 6);
  double r = 7;

  cout << "q0:      " << q0 << endl;
  cout << "q1:      " << q1 << endl;
  cout << "q2:      " << q2 << endl;
  cout << "r:       " << r << endl;
  cout << endl;
  cout << "-q0:     " << -q0 << endl;
  cout << "~q0:     " << ~q0 << endl;
  cout << endl;
  cout << "r * q0:  " << r*q0 << endl;
  cout << "r + q0:  " << r+q0 << endl;
  cout << "q0 / r:  " << q0/r << endl;
  cout << "q0 - r:  " << q0-r << endl;
  cout << endl;
  cout << "q0 + q1: " << q0+q1 << endl;
  cout << "q0 - q1: " << q0-q1 << endl;
  cout << "q0 * q1: " << q0*q1 << endl;
  cout << "q0 / q1: " << q0/q1 << endl;
  cout << endl;
  cout << "q0 * ~q0:     " << q0*~q0 << endl;
  cout << "q0 + q1*q2:   " << q0+q1*q2 << endl;
  cout << "(q0 + q1)*q2: " << (q0+q1)*q2 << endl;
  cout << "q0*q1*q2:     " << q0*q1*q2 << endl;
  cout << "(q0*q1)*q2:   " << (q0*q1)*q2 << endl;
  cout << "q0*(q1*q2):   " << q0*(q1*q2) << endl;
  cout << endl;
  cout << "||q0||:  " << sqrt(q0.normSquared()) << endl;
  cout << endl;
  cout << "q0*q1 - q1*q0: " << (q0*q1 - q1*q0) << endl;

  // Other base types
  Quaternion<int> q5(2), q6(3);
  cout << endl << q5*q6 << endl;
}
Output:
q0:      1 + 2i + 3j + 4k
q1:      2 + 3i + 4j + 5k
q2:      3 + 4i + 5j + 6k
r:       7

-q0:     -1 - 2i - 3j - 4k
~q0:     1 - 2i - 3j - 4k

r * q0:  7 + 14i + 21j + 28k
r + q0:  8 + 2i + 3j + 4k
q0 / r:  0.142857 + 0.285714i + 0.428571j + 0.571429k
q0 - r:  -6 + 2i + 3j + 4k

q0 + q1: 3 + 5i + 7j + 9k
q0 - q1: -1 - 1i - 1j - 1k
q0 * q1: -36 + 6i + 12j + 12k
q0 / q1: 0.740741 + 0i + 0.0740741j + 0.037037k

q0 * ~q0:     30 + 0i + 0j + 0k
q0 + q1*q2:   -55 + 18i + 27j + 30k
(q0 + q1)*q2: -100 + 24i + 42j + 42k
q0*q1*q2:     -264 - 114i - 132j - 198k
(q0*q1)*q2:   -264 - 114i - 132j - 198k
q0*(q1*q2):   -264 - 114i - 132j - 198k

||q0||:  5.47723

q0*q1 - q1*q0: 0 - 2i + 4j - 2k

6 + 0i + 0j + 0k

CLU

quat = cluster is make, minus, norm, conj, add, addr, mul, mulr,
                  equal, get_a, get_b, get_c, get_d, q_form 
    rep = struct[a,b,c,d: real]
    
    make = proc (a,b,c,d: real) returns (cvt)
        return (rep${a:a, b:b, c:c, d:d})
    end make
    
    minus = proc (q: cvt) returns (cvt)
        return (down(make(-q.a, -q.b, -q.c, -q.d)))
    end minus
    
    norm = proc (q: cvt) returns (real)
        return ((q.a**2.0 + q.b**2.0 + q.c**2.0 + q.d**2.0) ** 0.5)
    end norm
    
    conj = proc (q: cvt) returns (cvt)
        return (down(make(q.a, -q.b, -q.c, q.d)))
    end conj
    
    add = proc (q1, q2: cvt) returns (cvt)
        return (down(make(q1.a+q2.a, q1.b+q2.b, q1.c+q2.c, q1.d+q2.d)))
    end add
    
    addr = proc (q: cvt, r: real) returns (cvt)
        return (down(make(q.a+r, q.b+r, q.c+r, q.d+r)))
    end addr
    
    mul = proc (q1, q2: cvt) returns (cvt)
        a: real := q1.a*q2.a - q1.b*q2.b - q1.c*q2.c - q1.d*q2.d
        b: real := q1.a*q2.b + q1.b*q2.a + q1.c*q2.d - q1.d*q2.c
        c: real := q1.a*q2.c - q1.b*q2.d + q1.c*q2.a + q1.d*q2.b
        d: real := q1.a*q2.d + q1.b*q2.c - q1.c*q2.b + q1.d*q2.a
        return (down(make(a,b,c,d)))
    end mul
    
    mulr = proc (q: cvt, r: real) returns (cvt)
        return (down(make(q.a*r, q.b*r, q.c*r, q.d*r)))
    end mulr
    
    equal = proc (q1, q2: cvt) returns (bool)
        return (q1.a = q2.a & q1.b = q2.b & q1.c = q2.c & q1.d = q2.d)
    end equal
    
    get_a = proc (q: cvt) returns (real) return (q.a) end get_a 
    get_b = proc (q: cvt) returns (real) return (q.b) end get_b 
    get_c = proc (q: cvt) returns (real) return (q.c) end get_c 
    get_d = proc (q: cvt) returns (real) return (q.d) end get_d

    q_form = proc (q: cvt, a, b: int) returns (string)
        return ( f_form(q.a, a, b) || " + "
              || f_form(q.b, a, b) || "i + "
              || f_form(q.c, a, b) || "j + "
              || f_form(q.d, a, b) || "k" )
    end q_form
end quat
    
start_up = proc ()  
    po: stream := stream$primary_output()
    
    q0: quat := quat$make(1.0, 2.0, 3.0, 4.0)
    q1: quat := quat$make(2.0, 3.0, 4.0, 5.0)
    q2: quat := quat$make(3.0, 4.0, 5.0, 6.0)
    r:  real := 7.0
    
    stream$putl(po, "      q0 = " || quat$q_form(q0, 3, 3))
    stream$putl(po, "      q1 = " || quat$q_form(q1, 3, 3))
    stream$putl(po, "      q2 = " || quat$q_form(q2, 3, 3))
    stream$putl(po, "       r = " || f_form(r, 3, 3))
    stream$putl(po, "")
    
    stream$putl(po, "norm(q0) = " || f_form(quat$norm(q0), 3, 3))
    stream$putl(po, "     -q0 = " || quat$q_form(-q0, 3, 3))
    stream$putl(po, "conj(q0) = " || quat$q_form(quat$conj(q0), 3, 3))
    stream$putl(po, "  q0 + r = " || quat$q_form(quat$addr(q0, r), 3, 3))
    stream$putl(po, " q1 + q2 = " || quat$q_form(q1 + q2, 3, 3))
    stream$putl(po, "  q0 * r = " || quat$q_form(quat$mulr(q0, r), 3, 3))
    stream$putl(po, " q1 * q2 = " || quat$q_form(q1 * q2, 3, 3))
    stream$putl(po, " q2 * q1 = " || quat$q_form(q2 * q1, 3, 3))
    
    if q1*q2 ~= q2*q1 then stream$putl(po, "q1 * q2 ~= q2 * q1") end
end start_up
Output:
      q0 = 1.000 + 2.000i + 3.000j + 4.000k
      q1 = 2.000 + 3.000i + 4.000j + 5.000k
      q2 = 3.000 + 4.000i + 5.000j + 6.000k
       r = 7.000

norm(q0) = 5.477
     -q0 = -1.000 + -2.000i + -3.000j + -4.000k
conj(q0) = 1.000 + -2.000i + -3.000j + 4.000k
  q0 + r = 8.000 + 9.000i + 10.000j + 11.000k
 q1 + q2 = 5.000 + 7.000i + 9.000j + 11.000k
  q0 * r = 7.000 + 14.000i + 21.000j + 28.000k
 q1 * q2 = -56.000 + 16.000i + 24.000j + 26.000k
 q2 * q1 = -56.000 + 18.000i + 20.000j + 28.000k
q1 * q2 ~= q2 * q1

Common Lisp

(defclass quaternion () ((a :accessor q-a :initarg :a :type real)
                         (b :accessor q-b :initarg :b :type real)
                         (c :accessor q-c :initarg :c :type real)
                         (d :accessor q-d :initarg :d :type real))
  (:default-initargs :a 0 :b 0 :c 0 :d 0)) 

(defun make-q (&optional (a 0) (b 0) (c 0) (d 0))
  (make-instance 'quaternion :a a :b b :c c :d d))

(defgeneric sum (x y))

(defmethod sum ((x quaternion) (y quaternion))
  (make-q  (+ (q-a x) (q-a y))
           (+ (q-b x) (q-b y))
           (+ (q-c x) (q-c y))
           (+ (q-d x) (q-d y))))

(defmethod sum ((x quaternion) (y real))
  (make-q  (+ (q-a x) y) (q-b x) (q-c x) (q-d x)))

(defmethod sum ((x real) (y quaternion))
  (make-q  (+ (q-a y) x) (q-b y) (q-c y) (q-d y)))

(defgeneric sub (x y))

(defmethod sub ((x quaternion) (y quaternion))
  (make-q  (- (q-a x) (q-a y))
           (- (q-b x) (q-b y))
           (- (q-c x) (q-c y))
           (- (q-d x) (q-d y))))

(defmethod sub ((x quaternion) (y real))
  (make-q  (- (q-a x) y)
           (q-b x)
           (q-c x)
           (q-d x)))

(defmethod sub ((x real) (y quaternion))
  (make-q  (- (q-a y) x)
           (q-b y)
           (q-c y)
           (q-d y)))

(defgeneric mul (x y))

(defmethod mul ((x quaternion) (y real))
  (make-q  (* (q-a x) y)
           (* (q-b x) y)
           (* (q-c x) y)
           (* (q-d x) y)))

(defmethod mul ((x real) (y quaternion))
  (make-q  (* (q-a y) x)
           (* (q-b y) x)
           (* (q-c y) x)
           (* (q-d y) x)))

(defmethod mul ((x quaternion) (y quaternion))
  (make-q  (- (* (q-a x) (q-a y)) (* (q-b x) (q-b y)) (* (q-c x) (q-c y)) (* (q-d x) (q-d y)))
           (- (+ (* (q-a x) (q-b y)) (* (q-b x) (q-a y)) (* (q-c x) (q-d y))) (* (q-d x) (q-c y)))
           (- (+ (* (q-a x) (q-c y)) (* (q-c x) (q-a y)) (* (q-d x) (q-b y))) (* (q-b x) (q-d y)))
           (- (+ (* (q-a x) (q-d y)) (* (q-b x) (q-c y)) (* (q-d x) (q-a y))) (* (q-c x) (q-b y)))))

(defmethod norm ((x quaternion))
  (+ (sqrt (q-a x)) (sqrt (q-b x)) (sqrt (q-c x)) (sqrt (q-d x))))

(defmethod print-object ((x quaternion) stream)
  (format stream "~@f~@fi~@fj~@fk" (q-a x) (q-b x) (q-c x) (q-d x)))

(defvar q (make-q 0 1 0 0))
(defvar q1 (make-q 0 0 1 0))
(defvar q2 (make-q 0 0 0 1))
(defvar r 7)
(format t "q+q1+q2 = ~a~&" (reduce #'sum (list q q1 q2)))
(format t "r*(q+q1+q2) = ~a~&" (mul r (reduce #'sum (list q q1 q2))))
(format t "q*q1*q2 = ~a~&" (reduce #'mul (list q q1 q2)))
(format t "q-q1-q2 = ~a~&" (reduce #'sub (list q q1 q2)))
Output:
q+q1+q2 = +0.0+1.0i+1.0j+1.0k
r*(q+q1+q2) = +0.0+7.0i+7.0j+7.0k
q*q1*q2 = -1.0+0.0i+0.0j+0.0k
q-q1-q2 = +0.0+1.0i-1.0j-1.0k

Crystal

Translation of: Rust and Ruby
class Quaternion
  property a, b, c, d

  def initialize(@a : Int64, @b : Int64, @c : Int64, @d : Int64) end

  def norm; Math.sqrt(a**2 + b**2 + c**2 + d**2) end
  def conj; Quaternion.new(a, -b, -c, -d)        end
  def +(n)  Quaternion.new(a + n, b, c, d)       end
  def -(n)  Quaternion.new(a - n, b, c, d)       end
  def -()   Quaternion.new(-a, -b, -c, -d)       end
  def *(n)  Quaternion.new(a * n, b * n, c * n, d * n) end
  def ==(rhs : Quaternion) self.to_s == rhs.to_s end
  def +(rhs : Quaternion)
    Quaternion.new(a + rhs.a, b + rhs.b, c + rhs.c, d + rhs.d)
  end

  def -(rhs : Quaternion)
    Quaternion.new(a - rhs.a, b - rhs.b, c - rhs.c, d - rhs.d)
  end

  def *(rhs : Quaternion)
    Quaternion.new(
      a * rhs.a - b * rhs.b - c * rhs.c - d * rhs.d,
      a * rhs.b + b * rhs.a + c * rhs.d - d * rhs.c,
      a * rhs.c - b * rhs.d + c * rhs.a + d * rhs.b,
      a * rhs.d + b * rhs.c - c * rhs.b + d * rhs.a)
  end

  def to_s(io : IO) io << "(#{a} #{sgn(b)}i #{sgn(c)}j #{sgn(d)}k)\n" end
  private def sgn(n)  n.sign|1 == 1 ? "+ #{n}" : "- #{n.abs}"  end
end

struct Number
  def +(rhs : Quaternion)
    Quaternion.new(rhs.a + self, rhs.b, rhs.c, rhs.d)
  end

  def -(rhs : Quaternion)
    Quaternion.new(-rhs.a + self, -rhs.b, -rhs.c, -rhs.d)
  end

  def *(rhs : Quaternion)
    Quaternion.new(rhs.a * self, rhs.b * self, rhs.c * self, rhs.d * self)
  end
end

q0 = Quaternion.new(a: 1, b: 2, c: 3, d: 4)
q1 = Quaternion.new(2, 3, 4, 5)
q2 = Quaternion.new(3, 4, 5, 6)
r  = 7

puts "q0 = #{q0}"
puts "q1 = #{q1}"
puts "q2 = #{q2}"
puts "r  = #{r}"
puts
puts "normal of q0 = #{q0.norm}"
puts "-q0 = #{-q0}"
puts "conjugate of q0 = #{q0.conj}"
puts "q0 * (conjugate of q0) = #{q0 * q0.conj}"
puts "(conjugate of q0) * q0 = #{q0.conj * q0}"
puts
puts "r + q0 = #{r + q0}"
puts "q0 + r = #{q0 + r}"
puts
puts " q0 - r = #{q0 - r}"
puts "-q0 - r = #{-q0 - r}"
puts " r - q0 = #{r - q0}"
puts "-q0 + r = #{-q0 + r}"
puts
puts "r * q0 = #{r * q0}"
puts "q0 * r = #{q0 * r}"
puts
puts "q0 + q1 = #{q0 + q1}"
puts "q0 - q1 = #{q2 - q1}"
puts "q0 * q1 = #{q0 * q1}"
puts
puts " q0 + q1 * q2  = #{q0 + q1 * q2}"
puts "(q0 + q1) * q2 = #{(q0 + q1) * q2}"
puts
puts " q0 *  q1  * q2 = #{q0 * q1 * q2}"
puts "(q0 *  q1) * q2 = #{(q0 * q1) * q2}"
puts " q0 * (q1 * q2) = #{q0 * (q1 * q2)}"
puts
puts "q1 * q2 = #{q1 * q2}"
puts "q2 * q1 = #{q2 * q1}"
puts
puts "q1 * q2 != q2 * q1 => #{(q1 * q2) != (q2 * q1)}"
puts "q1 * q2 == q2 * q1 => #{(q1 * q2) == (q2 * q1)}"
Output:
q0 = (1 + 2i + 3j + 4k)
q1 = (2 + 3i + 4j + 5k)
q2 = (3 + 4i + 5j + 6k)
r  = 7

normal of q0 = 5.477225575051661
-q0 = (-1 - 2i - 3j - 4k)
conjugate of q0 = (1 - 2i - 3j - 4k)
q0 * (conjugate of q0) = (30 + 0i + 0j + 0k)
(conjugate of q0) * q0 = (30 + 0i + 0j + 0k)

r + q0 = (8 + 2i + 3j + 4k)
q0 + r = (8 + 2i + 3j + 4k)

 q0 - r = (-6 + 2i + 3j + 4k)
-q0 - r = (-8 - 2i - 3j - 4k)
 r - q0 = (6 - 2i - 3j - 4k)
-q0 + r = (6 - 2i - 3j - 4k)

r * q0 = (7 + 14i + 21j + 28k)
q0 * r = (7 + 14i + 21j + 28k)

q0 + q1 = (3 + 5i + 7j + 9k)
q0 - q1 = (1 + 1i + 1j + 1k)
q0 * q1 = (-36 + 6i + 12j + 12k)

 q0 + q1 * q2  = (-55 + 18i + 27j + 30k)
(q0 + q1) * q2 = (-100 + 24i + 42j + 42k)

 q0 *  q1  * q2 = (-264 - 114i - 132j - 198k)
(q0 *  q1) * q2 = (-264 - 114i - 132j - 198k)
 q0 * (q1 * q2) = (-264 - 114i - 132j - 198k)

q1 * q2 = (-56 + 16i + 24j + 26k)
q2 * q1 = (-56 + 18i + 20j + 28k)

q1 * q2 != q2 * q1 => true
q1 * q2 == q2 * q1 => false

D

import std.math, std.numeric, std.traits, std.conv, std.complex;


struct Quat(T) if (isFloatingPoint!T) {
    alias CT = Complex!T;

    union {
        struct { T re, i, j, k; } // Default init to NaN.
        struct { CT x, y; }
        struct { T[4] vector; }
    }

    string toString() const pure /*nothrow*/ @safe {
        return vector.text;
    }

    @property T norm2() const pure nothrow @safe @nogc { /// Norm squared.
        return re ^^ 2 + i ^^ 2 + j ^^ 2 + k ^^ 2;
    }

    @property T abs() const pure nothrow @safe @nogc { /// Norm.
        return sqrt(norm2);
    }

    @property T arg() const pure nothrow @safe @nogc { /// Theta.
        return acos(re / abs); // this may be incorrect...
    }

    @property Quat!T conj() const pure nothrow @safe @nogc { /// Conjugate.
        return Quat!T(re, -i, -j, -k);
    }

    @property Quat!T recip() const pure nothrow @safe @nogc {  /// Reciprocal.
        return Quat!T(re / norm2, -i / norm2, -j / norm2, -k / norm2);
    }

    @property Quat!T pureim() const pure nothrow @safe @nogc { /// Pure imagery.
        return Quat!T(0, i, j, k);
    }

    @property Quat!T versor() const pure nothrow @safe @nogc { /// Unit versor.
        return this / abs;
    }

    /// Unit versor of imagery part.
    @property Quat!T iversor() const pure nothrow @safe @nogc {
        return pureim / pureim.abs;
    }

    /// Assignment.
    Quat!T opAssign(U : T)(Quat!U z) pure nothrow @safe @nogc {
        x = z.x;  y = z.y;
        return this;
    }

    Quat!T opAssign(U : T)(Complex!U c) pure nothrow @safe @nogc {
        x = c;  y = 0;
        return this;
    }

    Quat!T opAssign(U : T)(U r) pure nothrow @safe @nogc
    if (isNumeric!U) {
        re = r; i = 0; y = 0;
        return this;
    }

    /// Test for equal, not ordered so no opCmp.
    bool opEquals(U : T)(Quat!U z) const pure nothrow @safe @nogc {
        return re == z.re && i == z.i && j == z.j && k == z.k;
    }

    bool opEquals(U : T)(Complex!U c) const pure nothrow @safe @nogc {
        return re == c.re && i == c.im && j == 0 && k == 0;
    }

    bool opEquals(U : T)(U r) const pure nothrow @safe @nogc
    if (isNumeric!U) {
        return re == r && i == 0 && j == 0 && k == 0;
    }

    /// Unary op.
    Quat!T opUnary(string op)() const pure nothrow @safe @nogc
    if (op == "+") {
        return this;
    }

    Quat!T opUnary(string op)() const pure nothrow @safe @nogc
    if (op == "-") {
        return Quat!T(-re, -i, -j, -k);
    }

    /// Binary op, Quaternion on left of op.
    Quat!(CommonType!(T,U)) opBinary(string op, U)(Quat!U z)
    const pure nothrow @safe @nogc {
        alias typeof(return) C;

        static if (op == "+" ) {
            return C(re + z.re, i + z.i, j + z.j, k + z.k);
        } else static if (op == "-") {
            return C(re - z.re, i - z.i, j - z.j, k - z.k);
        } else static if (op == "*") {
            return C(re * z.re - i * z.i  - j * z.j  - k * z.k,
                     re * z.i  + i * z.re + j * z.k  - k * z.j,
                     re * z.j  - i * z.k  + j * z.re + k * z.i,
                     re * z.k  + i * z.j  - j * z.i  + k * z.re);
        } else static if (op == "/") {
            return this * z.recip;
        }
    }

    /// Extend complex to quaternion.
    Quat!(CommonType!(T,U)) opBinary(string op, U)(Complex!U c)
    const pure nothrow @safe @nogc {
        return opBinary!op(typeof(return)(c.re, c.im, 0, 0));
    }

    /// For scalar.
    Quat!(CommonType!(T,U)) opBinary(string op, U)(U r)
    const pure nothrow @safe @nogc
    if (isNumeric!U) {
        alias typeof(return) C;

        static if (op == "+" ) {
            return C(re + r, i, j, k);
        } else static if (op == "-") {
            return C(re - r, i, j, k);
        } else static if (op == "*") {
            return C(re * r, i * r, j * r, k * r);
        } else static if (op == "/") {
            return C(re / r, i / r, j / r, k / r);
        } else static if (op == "^^") {
            return pow(r);
        }
    }

    /// Power function.
    Quat!(CommonType!(T,U)) pow(U)(U r)
    const pure nothrow @safe @nogc
    if (isNumeric!U) {
        return (abs^^r) * exp(r * iversor * arg);
    }

    /// Handle binary op if Quaternion on right of op and left is
    /// not quaternion.
    Quat!(CommonType!(T,U)) opBinaryRight(string op, U)(Complex!U c)
    const pure nothrow @safe @nogc {
        alias typeof(return) C;
        auto w = C(c.re, c.im, 0, 0);
        return w.opBinary!(op)(this);
    }

    Quat!(CommonType!(T,U)) opBinaryRight(string op, U)(U r)
    const pure nothrow @safe @nogc
    if (isNumeric!U) {
        alias typeof(return) C;

        static if (op == "+" || op == "*") {
            return opBinary!op(r);
        } else static if (op == "-") {
            return C(r - re , -i, -j, -k);
        } else static if (op == "/") {
            auto w = C(re, i, j, k);
            return w.recip * r;
        }
    }
}


HT exp(HT)(HT z) pure nothrow @safe @nogc
if (is(HT T == Quat!T)) {
    immutable inorm = z.pureim.abs;
    return std.math.exp(z.re) * (cos(inorm) + z.iversor * sin(inorm));
}

HT log(HT)(HT z) pure nothrow @safe @nogc
if (is(HT T == Quat!T)) {
    return std.math.log(z.abs) + z.iversor * acos(z.re / z.abs);
}


void main() @safe { // Demo code.
    import std.stdio;

    alias QR = Quat!real;
    enum real r = 7.0;

    immutable QR q  = QR(2, 3, 4, 5),
                 q1 = QR(2, 3, 4, 5),
                 q2 = QR(3, 4, 5, 6);

    writeln("1.             q - norm: ", q.abs);
    writeln("2.         q - negative: ", -q);
    writeln("3.        q - conjugate: ", q.conj);
    writeln("4.                r + q: ", r + q);
    writeln("                  q + r: ", q + r);
    writeln("5.              q1 + q2: ", q1 + q2);
    writeln("6.                r * q: ", r * q);
    writeln("                  q * r: ", q * r);
    writeln("7.              q1 * q2: ", q1 * q2);
    writeln("                q2 * q1: ", q2 * q1);
    writeln("8.  q1 * q2 != q2 * Q1 ? ", q1 * q2 != q2 * q1);

    immutable QR i = QR(0, 1, 0, 0),
                 j = QR(0, 0, 1, 0),
                 k = QR(0, 0, 0, 1);
    writeln("9.1               i * i: ", i * i);
    writeln("                  J * j: ", j * j);
    writeln("                  k * k: ", k * k);
    writeln("              i * j * k: ", i * j * k);
    writeln("9.2             q1 / q2: ", q1 / q2);
    writeln("9.3        q1 / q2 * q2: ", q1 / q2 * q2);
    writeln("           q2 * q1 / q2: ", q2 * q1 / q2);
    writeln("9.4         exp(pi * i): ", exp(PI * i));
    writeln("            exp(pi * j): ", exp(PI * j));
    writeln("            exp(pi * k): ", exp(PI * k));
    writeln("                 exp(q): ", exp(q));
    writeln("                 log(q): ", log(q));
    writeln("            exp(log(q)): ", exp(log(q)));
    writeln("            log(exp(q)): ", log(exp(q)));
    immutable s = q.exp.log;
    writeln("9.5 let s = log(exp(q)): ", s);
    writeln("                 exp(s): ", exp(s));
    writeln("                 log(s): ", log(s));
    writeln("            exp(log(s)): ", exp(log(s)));
    writeln("            log(exp(s)): ", log(exp(s)));
}
Output:
1.             q - norm: 7.34847
2.         q - negative: [-2, -3, -4, -5]
3.        q - conjugate: [2, -3, -4, -5]
4.                r + q: [9, 3, 4, 5]
                  q + r: [9, 3, 4, 5]
5.              q1 + q2: [5, 7, 9, 11]
6.                r * q: [14, 21, 28, 35]
                  q * r: [14, 21, 28, 35]
7.              q1 * q2: [-56, 16, 24, 26]
                q2 * q1: [-56, 18, 20, 28]
8.  q1 * q2 != q2 * Q1 ? true
9.1               i * i: [-1, 0, 0, 0]
                  J * j: [-1, 0, 0, 0]
                  k * k: [-1, 0, 0, 0]
              i * j * k: [-1, 0, 0, 0]
9.2             q1 / q2: [0.790698, 0.0232558, -1.35525e-20, 0.0465116]
9.3        q1 / q2 * q2: [2, 3, 4, 5]
           q2 * q1 / q2: [2, 3.46512, 3.90698, 4.76744]
9.4         exp(pi * i): [-1, -5.42101e-20, -0, -0]
            exp(pi * j): [-1, -0, -5.42101e-20, -0]
            exp(pi * k): [-1, -0, -0, -5.42101e-20]
                 exp(q): [5.21186, 2.22222, 2.96296, 3.7037]
                 log(q): [1.99449, 0.549487, 0.732649, 0.915812]
            exp(log(q)): [2, 3, 4, 5]
            log(exp(q)): [2, 0.33427, 0.445694, 0.557117]
9.5 let s = log(exp(q)): [2, 0.33427, 0.445694, 0.557117]
                 exp(s): [5.21186, 2.22222, 2.96296, 3.7037]
                 log(s): [0.765279, 0.159215, 0.212286, 0.265358]
            exp(log(s)): [2, 0.33427, 0.445694, 0.557117]
            log(exp(s)): [2, 0.33427, 0.445694, 0.557117]


Dart

Translation of: Kotlin
import 'dart:math' as math;

class Quaternion {
  final double a, b, c, d;

  Quaternion(this.a, this.b, this.c, this.d);

  Quaternion operator +(Object other) {
    if (other is Quaternion) {
      return Quaternion(a + other.a, b + other.b, c + other.c, d + other.d);
    } else if (other is double) {
      return Quaternion(a + other, b, c, d);
    }
    throw ArgumentError('Invalid type for addition: ${other.runtimeType}');
  }

  Quaternion operator *(Object other) {
    if (other is Quaternion) {
      return Quaternion(
        a * other.a - b * other.b - c * other.c - d * other.d,
        a * other.b + b * other.a + c * other.d - d * other.c,
        a * other.c - b * other.d + c * other.a + d * other.b,
        a * other.d + b * other.c - c * other.b + d * other.a,
      );
    } else if (other is double) {
      return Quaternion(a * other, b * other, c * other, d * other);
    }
    throw ArgumentError('Invalid type for multiplication: ${other.runtimeType}');
  }

  Quaternion operator -() => Quaternion(-a, -b, -c, -d);

  Quaternion conj() => Quaternion(a, -b, -c, -d);

  double norm() => math.sqrt(a * a + b * b + c * c + d * d);

  @override
  String toString() => '($a, $b, $c, $d)';
}

void main() {
  var q  = Quaternion(1.0, 2.0, 3.0, 4.0);
  var q1 = Quaternion(2.0, 3.0, 4.0, 5.0);
  var q2 = Quaternion(3.0, 4.0, 5.0, 6.0);
  var r  = 7.0;
  print("q  = $q");
  print("q1 = $q1");
  print("q2 = $q2");
  print("r  = $r\n");
  print("norm(q) = ${q.norm().toStringAsFixed(6)}");
  print("-q      = ${-q}");
  print("conj(q) = ${q.conj()}\n");
  print("r  + q  = ${q + r}");
  print("q  + r  = ${q + r}");
  print("q1 + q2 = ${q1 + q2}\n");
  print("r  * q  = ${q * r}");
  print("q  * r  = ${q * r}");
  var q3 = q1 * q2;
  var q4 = q2 * q1;
  print("q1 * q2 = $q3");
  print("q2 * q1 = $q4\n");
  print("q1 * q2 != q2 * q1 = ${q3 != q4}");
}
Output:
q  = (1.0, 2.0, 3.0, 4.0)
q1 = (2.0, 3.0, 4.0, 5.0)
q2 = (3.0, 4.0, 5.0, 6.0)
r  = 7.0

norm(q) = 5.477226
-q      = (-1.0, -2.0, -3.0, -4.0)
conj(q) = (1.0, -2.0, -3.0, -4.0)

r  + q  = (8.0, 2.0, 3.0, 4.0)
q  + r  = (8.0, 2.0, 3.0, 4.0)
q1 + q2 = (5.0, 7.0, 9.0, 11.0)

r  * q  = (7.0, 14.0, 21.0, 28.0)
q  * r  = (7.0, 14.0, 21.0, 28.0)
q1 * q2 = (-56.0, 16.0, 24.0, 26.0)
q2 * q1 = (-56.0, 18.0, 20.0, 28.0)

q1 * q2 != q2 * q1 = true


Delphi

unit Quaternions;

interface

type

  TQuaternion = record
    A, B, C, D: double;

    function  Init          (aA, aB, aC, aD : double): TQuaternion;
    function  Norm          : double;
    function  Conjugate     : TQuaternion;
    function  ToString      : string;

    class operator Negative (Left : TQuaternion): TQuaternion;
    class operator Positive (Left : TQuaternion): TQuaternion;
    class operator Add      (Left, Right : TQuaternion): TQuaternion;
    class operator Add      (Left : TQuaternion; Right : double): TQuaternion; overload;
    class operator Add      (Left : double; Right : TQuaternion): TQuaternion; overload;
    class operator Subtract (Left, Right : TQuaternion): TQuaternion;
    class operator Multiply (Left, Right : TQuaternion): TQuaternion;
    class operator Multiply (Left : TQuaternion; Right : double): TQuaternion; overload;
    class operator Multiply (Left : double; Right : TQuaternion): TQuaternion; overload;
  end;

implementation

uses
  SysUtils;

{ TQuaternion }

function TQuaternion.Init(aA, aB, aC, aD: double): TQuaternion;
begin
  A := aA;
  B := aB;
  C := aC;
  D := aD;

  result := Self;
end;

function TQuaternion.Norm: double;
begin
  result := sqrt(sqr(A) + sqr(B) + sqr(C) + sqr(D));
end;

function TQuaternion.Conjugate: TQuaternion;
begin
  result.B := -B;
  result.C := -C;
  result.D := -D;
end;

class operator TQuaternion.Negative(Left: TQuaternion): TQuaternion;
begin
  result.A := -Left.A;
  result.B := -Left.B;
  result.C := -Left.C;
  result.D := -Left.D;
end;

class operator TQuaternion.Positive(Left: TQuaternion): TQuaternion;
begin
  result := Left;
end;

class operator TQuaternion.Add(Left, Right: TQuaternion): TQuaternion;
begin
  result.A := Left.A + Right.A;
  result.B := Left.B + Right.B;
  result.C := Left.C + Right.C;
  result.D := Left.D + Right.D;
end;

class operator TQuaternion.Add(Left: TQuaternion; Right: double): TQuaternion;
begin
  result.A := Left.A + Right;
  result.B := Left.B;
  result.C := Left.C;
  result.D := Left.D;
end;

class operator TQuaternion.Add(Left: double; Right: TQuaternion): TQuaternion;
begin
  result.A := Left + Right.A;
  result.B := Right.B;
  result.C := Right.C;
  result.D := Right.D;
end;

class operator TQuaternion.Subtract(Left, Right: TQuaternion): TQuaternion;
begin
  result.A := Left.A - Right.A;
  result.B := Left.B - Right.B;
  result.C := Left.C - Right.C;
  result.D := Left.D - Right.D;
end;

class operator TQuaternion.Multiply(Left, Right: TQuaternion): TQuaternion;
begin
  result.A := Left.A * Right.A - Left.B * Right.B - Left.C * Right.C - Left.D * Right.D;
  result.B := Left.A * Right.B + Left.B * Right.A + Left.C * Right.D - Left.D * Right.C;
  result.C := Left.A * Right.C - Left.B * Right.D + Left.C * Right.A + Left.D * Right.B;
  result.D := Left.A * Right.D + Left.B * Right.C - Left.C * Right.B + Left.D * Right.A;
end;

class operator TQuaternion.Multiply(Left: double; Right: TQuaternion): TQuaternion;
begin
  result.A := Left * Right.A;
  result.B := Left * Right.B;
  result.C := Left * Right.C;
  result.D := Left * Right.D;
end;

class operator TQuaternion.Multiply(Left: TQuaternion; Right: double): TQuaternion;
begin
  result.A := Left.A * Right;
  result.B := Left.B * Right;
  result.C := Left.C * Right;
  result.D := Left.D * Right;
end;

function TQuaternion.ToString: string;
begin
  result := Format('%f + %fi + %fj + %fk', [A, B, C, D]);
end;

end.

Test program

program QuaternionTest;

{$APPTYPE CONSOLE}

uses
  Quaternions in 'Quaternions.pas';

var
  r : double;
  q, q1, q2 : TQuaternion;
begin
  r := 7;
  q  := q .Init(1, 2, 3, 4);
  q1 := q1.Init(2, 3, 4, 5);
  q2 := q2.Init(3, 4, 5, 6);

  writeln('q  = ', q.ToString);
  writeln('q1 = ', q1.ToString);
  writeln('q2 = ', q2.ToString);
  writeln('r  = ', r);
  writeln('Norm(q ) = ', q.Norm);
  writeln('Norm(q1) = ', q1.Norm);
  writeln('Norm(q2) = ', q2.Norm);
  writeln('-q = ', (-q).ToString);
  writeln('Conjugate q = ', q.Conjugate.ToString);
  writeln('q1 + q2 = ', (q1 + q2).ToString);
  writeln('q2 + q1 = ', (q2 + q1).ToString);
  writeln('q * r   = ', (q * r).ToString);
  writeln('r * q   = ', (r * q).ToString);
  writeln('q1 * q2 = ', (q1 * q2).ToString);
  writeln('q2 * q1 = ', (q2 * q1).ToString);
end.
Output:
q  = 1.00 + 2.00i + 3.00j + 4.00k
q1 = 2.00 + 3.00i + 4.00j + 5.00k
q2 = 3.00 + 4.00i + 5.00j + 6.00k
r  =  7.00000000000000E+0000
Norm(q ) =  5.47722557505166E+0000
Norm(q1) =  7.34846922834953E+0000
Norm(q2) =  9.27361849549570E+0000
-q = -1.00 + -2.00i + -3.00j + -4.00k
Conjugate q = -1.00 + -2.00i + -3.00j + -4.00k
q1 + q2 = 5.00 + 7.00i + 9.00j + 11.00k
q2 + q1 = 5.00 + 7.00i + 9.00j + 11.00k
q * r   = 7.00 + 14.00i + 21.00j + 28.00k
r * q   = 7.00 + 14.00i + 21.00j + 28.00k
q1 * q2 = -56.00 + 16.00i + 24.00j + 26.00k
q2 * q1 = -56.00 + 18.00i + 20.00j + 28.00k

--DavidIzadaR 20:33, 7 August 2011 (UTC)

E

interface Quaternion guards QS {}
def makeQuaternion(a, b, c, d) {
    return def quaternion implements QS {
    
        to __printOn(out) {
            out.print("(", a, " + ", b, "i + ")
            out.print(c, "j + ", d, "k)")
        }

        # Task requirement 1
        to norm() {
            return (a**2 + b**2 + c**2 + d**2).sqrt()
        }

        # Task requirement 2
        to negate() {
            return makeQuaternion(-a, -b, -c, -d)
        }
        
        # Task requirement 3
        to conjugate() {
            return makeQuaternion(a, -b, -c, -d)
        }

        # Task requirement 4, 5
        # This implements q + r; r + q is deliberately prohibited by E
        to add(other :any[Quaternion, int, float64]) {
            switch (other) {
              match q :Quaternion {
                return makeQuaternion(
                    a+q.a(), b+q.b(), c+q.c(), d+q.d())
              }
              match real {
                return makeQuaternion(a+real, b, c, d)
              }
            }
        }

        # Task requirement 6, 7
        # This implements q * r; r * q is deliberately prohibited by E
        to multiply(other :any[Quaternion, int, float64]) {
            switch (other) {
                match q :Quaternion {
                    return makeQuaternion(
                        a*q.a() - b*q.b() - c*q.c() - d*q.d(),
                        a*q.b() + b*q.a() + c*q.d() - d*q.c(),
                        a*q.c() - b*q.d() + c*q.a() + d*q.b(),
                        a*q.d() + b*q.c() - c*q.b() + d*q.a())
                    }
                match real {
                    return makeQuaternion(real*a, real*b, real*c, real*d)
                }
            }
        }
        
        to a() { return a }
        to b() { return b }
        to c() { return c }
        to d() { return d }
    }
}
? def q1 := makeQuaternion(2,3,4,5)
# value: (2 + 3i + 4j + 5k)

? def q2 := makeQuaternion(3,4,5,6)
# value: (3 + 4i + 5j + 6k)

? q1+q2
# value: (5 + 7i + 9j + 11k)

? q1*q2
# value: (-56 + 16i + 24j + 26k)

? q2*q1
# value: (-56 + 18i + 20j + 28k)

? q1+(-2)
# value: (0 + 3i + 4j + 5k)

EasyLang

func qnorm q[] .
   for i to 4
      s += q[i] * q[i]
   .
   return sqrt s
.
func[] qneg q[] .
   for i to 4
      q[i] = -q[i]
   .
   return q[]
.
func[] qconj q[] .
   for i = 2 to 4
      q[i] = -q[i]
   .
   return q[]
.
func[] qaddreal q[] r .
   q[1] += r
   return q[]
.
func[] qadd q[] q2[] .
   for i to 4
      q[i] += q2[i]
   .
   return q[]
.
func[] qmulreal q[] r .
   for i to 4
      q[i] *= r
   .
   return q[]
.
func[] qmul q1[] q2[] .
   res[] &= q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3] - q1[4] * q2[4]
   res[] &= q1[1] * q2[2] + q1[2] * q2[1] + q1[3] * q2[4] - q1[4] * q2[3]
   res[] &= q1[1] * q2[3] - q1[2] * q2[4] + q1[3] * q2[1] + q1[4] * q2[2]
   res[] &= q1[1] * q2[4] + q1[2] * q2[3] - q1[3] * q2[2] + q1[4] * q2[1]
   return res[]
.
q[] = [ 1 2 3 4 ]
q1[] = [ 2 3 4 5 ]
q2[] = [ 3 4 5 6 ]
r = 7
# 
print "q = " & q[]
print "q1 = " & q1[]
print "q2 = " & q2[]
print "r = " & r
print "norm(q) = " & qnorm q[]
print "neg(q) = " & qneg q[]
print "conjugate(q) = " & qconj q[]
print "q+r = " & qaddreal q[] r
print "q1+q2 = " & qadd q1[] q2[]
print "qr = " & qmulreal q[] r
print "q1q2 = " & qmul q1[] q2[]
print "q2q1 = " & qmul q2[] q1[]
if q1[] <> q2[]
   print "q1 != q2"
.
Output:
q = [ 1 2 3 4 ]
q1 = [ 2 3 4 5 ]
q2 = [ 3 4 5 6 ]
r = 7
norm(q) = 5.48
neg(q) = [ -1 -2 -3 -4 ]
conjugate(q) = [ 1 -2 -3 -4 ]
q+r = [ 8 2 3 4 ]
q1+q2 = [ 5 7 9 11 ]
qr = [ 7 14 21 28 ]
q1q2 = [ -56 16 24 26 ]
q2q1 = [ -56 18 20 28 ]
q1 != q2

Eero

#import <Foundation/Foundation.h>

interface Quaternion : Number
  // Properties -- note that this is an immutable class.
  double real, i, j, k {readonly}
end

implementation Quaternion

  initWithReal: double, i: double, j: double, k: double, return instancetype
    self = super.init
    if self
      _real = real; _i = i; _j = j; _k = k
    return self

  +new: double real, ..., return instancetype
    va_list args
    va_start(args, real)
    object := Quaternion.alloc.initWithReal: real, 
                                          i: va_arg(args, double),
                                          j: va_arg(args, double),
                                          k: va_arg(args, double)
    va_end(args)
    return object

  descriptionWithLocale: id, return String = String.stringWithFormat:
      '(%.1f, %.1f, %.1f, %.1f)', self.real, self.i, self.j, self.k

  norm, return double =
      sqrt(self.real * self.real +
           self.i * self.i + self.j * self.j + self.k * self.k)

  negative, return Quaternion =
      Quaternion.new: -self.real, -self.i, -self.j, -self.k

  conjugate, return Quaternion =
      Quaternion.new: self.real, -self.i, -self.j, -self.k

  // Overload "+" operator (left operand is Quaternion)
  plus: Number operand, return Quaternion
    real := self.real, i = self.i, j = self.j, k = self.k
    if operand.isKindOfClass: Quaternion.class
      q := (Quaternion)operand
      real += q.real; i += q.i; j += q.j; k += q.k
    else
      real += (double)operand
    return Quaternion.new: real, i, j, k 

  // Overload "*" operator (left operand is Quaternion)
  multipliedBy: Number operand, return Quaternion
    real := self.real, i = self.i, j = self.j, k = self.k
    if operand.isKindOfClass: Quaternion.class
      q := (Quaternion)operand
      real = self.real * q.real - self.i* q.i - self.j * q.j - self.k * q.k
      i = self.real * q.i + self.i * q.real + self.j * q.k - self.k * q.j
      j = self.real * q.j - self.i * q.k + self.j * q.real + self.k * q.i
      k = self.real * q.k + self.i * q.j - self.j * q.i + self.k * q.real
    else
      real *= (double)operand
      i *= (double)operand; j *= (double)operand; k *= (double)operand
    return Quaternion.new: real, i, j, k

end

implementation Number (QuaternionOperators)

  // Overload "+" operator (left operand is Number)
  plus: Quaternion operand, return Quaternion
    real := (double)self + operand.real
    return Quaternion.new: real, operand.i, operand.j, operand.k

  // Overload "*" operator (left operand is Number)
  multipliedBy: Quaternion operand, return Quaternion
    r := (double)self
    return Quaternion.new: r * operand.real, r * operand.i,
                           r * operand.j, r * operand.k

end

int main()
  autoreleasepool

    q  := Quaternion.new: 1.0, 2.0, 3.0, 4.0
    q1 := Quaternion.new: 2.0, 3.0, 4.0, 5.0
    q2 := Quaternion.new: 3.0, 4.0, 5.0, 6.0

    Log( 'q  = %@', q )
    Log( 'q1 = %@', q1 )
    Log( 'q2 = %@\n\n', q2 )

    Log( 'q norm = %.3f',  q.norm )
    Log( 'q negative = %@',  q.negative )
    Log( 'q conjugate = %@',  q.conjugate )
    Log( '7 + q = %@', 7.0 + q )
    Log( 'q + 7 = %@', q + 7.0 )
    Log( 'q1 + q2 = %@',  q1 + q2 )
    Log( '7 * q = %@', 7 * q)
    Log( 'q * 7 = %@', q * 7.0 )
    Log( 'q1 * q2 = %@',  q1 * q2 )
    Log( 'q2 * q1 = %@',  q2 * q1 )

  return 0
Output:
2013-09-04 16:40:29.818 a.out[2170:507] q  = (1.0, 2.0, 3.0, 4.0)
2013-09-04 16:40:29.819 a.out[2170:507] q1 = (2.0, 3.0, 4.0, 5.0)
2013-09-04 16:40:29.820 a.out[2170:507] q2 = (3.0, 4.0, 5.0, 6.0)

2013-09-04 16:40:29.820 a.out[2170:507] q norm = 5.477
2013-09-04 16:40:29.820 a.out[2170:507] q negative = (-1.0, -2.0, -3.0, -4.0)
2013-09-04 16:40:29.820 a.out[2170:507] q conjugate = (1.0, -2.0, -3.0, -4.0)
2013-09-04 16:40:29.821 a.out[2170:507] 7 + q = (8.0, 2.0, 3.0, 4.0)
2013-09-04 16:40:29.821 a.out[2170:507] q + 7 = (8.0, 2.0, 3.0, 4.0)
2013-09-04 16:40:29.821 a.out[2170:507] q1 + q2 = (5.0, 7.0, 9.0, 11.0)
2013-09-04 16:40:29.821 a.out[2170:507] 7 * q = (7.0, 14.0, 21.0, 28.0)
2013-09-04 16:40:29.821 a.out[2170:507] q * 7 = (7.0, 14.0, 21.0, 28.0)
2013-09-04 16:40:29.822 a.out[2170:507] q1 * q2 = (-56.0, 16.0, 24.0, 26.0)
2013-09-04 16:40:29.822 a.out[2170:507] q2 * q1 = (-56.0, 18.0, 20.0, 28.0)

Elena

Translation of: C#

ELENA 6.x :

import system'math;
import extensions;
import extensions'text;
 
struct Quaternion
{
    real A : rprop;
    real B : rprop;
    real C : rprop;
    real D : rprop;
 
    constructor new(a, b, c, d)
        <= new(cast real(a), cast real(b), cast real(c), cast real(d));
 
    constructor new(real a, real b, real c, real d)
    {
        A := a;
        B := b;
        C := c;
        D := d
    }
 
    constructor(real r)
    {
        A := r;
        B := 0.0r;
        C := 0.0r;
        D := 0.0r
    }
 
    real Norm = (A*A + B*B + C*C + D*D).sqrt();
 
    Quaternion Negative = Quaternion.new(A.Negative,B.Negative,C.Negative,D.Negative);
 
    Quaternion Conjugate = Quaternion.new(A,B.Negative,C.Negative,D.Negative);
 
    Quaternion add(Quaternion q)
        = Quaternion.new(A + q.A, B + q.B, C + q.C, D + q.D);
 
    Quaternion multiply(Quaternion q)
        = Quaternion.new(
            A * q.A - B * q.B - C * q.C - D * q.D,
            A * q.B + B * q.A + C * q.D - D * q.C,
            A * q.C - B * q.D + C * q.A + D * q.B,
            A * q.D + B * q.C - C * q.B + D * q.A);
 
    Quaternion add(real r)
        <= add(Quaternion.new(r,0,0,0));
 
    Quaternion multiply(real r)
        <= multiply(Quaternion.new(r,0,0,0));
 
    bool equal(Quaternion q)
        = (A == q.A) && (B == q.B) && (C == q.C) && (D == q.D);
 
    string toPrintable()
        = new StringWriter().printFormatted("Q({0}, {1}, {2}, {3})",A,B,C,D);
}
 
public program()
{
    auto q := Quaternion.new(1,2,3,4);
    auto q1 := Quaternion.new(2,3,4,5);
    auto q2 := Quaternion.new(3,4,5,6);
    real r := 7;
 
    console.printLine("q = ", q);
    console.printLine("q1 = ", q1);
    console.printLine("q2 = ", q2);
    console.printLine("r = ", r);
 
    console.printLine("q.Norm() = ", q.Norm);
    console.printLine("q1.Norm() = ", q1.Norm);
    console.printLine("q2.Norm() = ", q2.Norm);
 
    console.printLine("-q = ", q.Negative);
    console.printLine("q.Conjugate() = ", q.Conjugate);
 
    console.printLine("q + r = ", q + r);
    console.printLine("q1 + q2 = ", q1 + q2);
    console.printLine("q2 + q1 = ", q2 + q1);
 
    console.printLine("q * r = ", q * r);
    console.printLine("q1 * q2 = ", q1 * q2);
    console.printLine("q2 * q1 = ", q2 * q1);
 
    console.printLineFormatted("q1*q2 {0} q2*q1", ((q1 * q2) == (q2 * q1)).iif("==","!="))
}
Output:
q = Q(1.0, 2.0, 3.0, 4.0)
q1 = Q(2.0, 3.0, 4.0, 5.0)
q2 = Q(3.0, 4.0, 5.0, 6.0)
r = 7.0
q.Norm() = 5.477225575052
q1.Norm() = 7.34846922835
q2.Norm() = 9.273618495496
-q = Q(-1.0, -2.0, -3.0, -4.0)
q.Conjugate() = Q(1.0, -2.0, -3.0, -4.0)
q + r = Q(8.0, 2.0, 3.0, 4.0)
q1 + q2 = Q(5.0, 7.0, 9.0, 11.0)
q2 + q1 = Q(5.0, 7.0, 9.0, 11.0)
q * r = Q(7.0, 14.0, 21.0, 28.0)
q1 * q2 = Q(-56.0, 16.0, 24.0, 26.0)
q2 * q1 = Q(-56.0, 18.0, 20.0, 28.0)
q1*q2 != q2*q1

ERRE

PROGRAM QUATERNION

!$DOUBLE

TYPE QUATERNION=(A,B,C,D)

DIM Q:QUATERNION,Q1:QUATERNION,Q2:QUATERNION


DIM R:QUATERNION,S:QUATERNION,T:QUATERNION

PROCEDURE NORM(T.->NORM)
   NORM=SQR(T.A*T.A+T.B*T.B+T.C*T.C+T.D*T.D)
END PROCEDURE

PROCEDURE NEGATIVE(T.->T.)
    T.A=-T.A
    T.B=-T.B
    T.C=-T.C
    T.D=-T.D
END PROCEDURE

PROCEDURE CONJUGATE(T.->T.)
    T.A=T.A
    T.B=-T.B
    T.C=-T.C
    T.D=-T.D
END PROCEDURE

PROCEDURE ADD_REAL(T.,REAL->T.)
    T.A=T.A+REAL
    T.B=T.B
    T.C=T.C
    T.D=T.D
END PROCEDURE

PROCEDURE ADD(T.,S.->T.)
    T.A=T.A+S.A
    T.B=T.B+S.B
    T.C=T.C+S.C
    T.D=T.D+S.D
END PROCEDURE

PROCEDURE MULT_REAL(T.,REAL->T.)
    T.A=T.A*REAL
    T.B=T.B*REAL
    T.C=T.C*REAL
    T.D=T.D*REAL
END PROCEDURE

PROCEDURE MULT(T.,S.->R.)
    R.A=T.A*S.A-T.B*S.B-T.C*S.C-T.D*S.D
    R.B=T.A*S.B+T.B*S.A+T.C*S.D-T.D*S.C
    R.C=T.A*S.C-T.B*S.D+T.C*S.A+T.D*S.B
    R.D=T.A*S.D+T.B*S.C-T.C*S.B+T.D*S.A
END PROCEDURE

PROCEDURE PRINTQ(T.)
    PRINT("(";T.A;",";T.B;",";T.C;",";T.D;")")
END PROCEDURE

BEGIN
    Q.A=1  Q.B=2  Q.C=3  Q.D=4
    Q1.A=2 Q1.B=3 Q1.C=4 Q1.D=5
    Q2.A=3 Q2.B=4 Q2.C=5 Q2.D=6
    REAL=7

    NORM(Q.->NORM)
    PRINT("Norm(q)=";NORM)

    NEGATIVE(Q.->T.)
    PRINT("Negative(q) =";)
    PRINTQ(T.)

    CONJUGATE(Q.->T.)
    PRINT("Conjugate(q) =";)
    PRINTQ(T.)

    ADD_REAL(Q.,REAL->T.)
    PRINT("q + real =";)
    PRINTQ(T.)

! addition is commutative
    ADD(Q1.,Q2.->T.)
    PRINT("q1 + q2 =";)
    PRINTQ(T.)

    ADD(Q2.,Q1.->T.)
    PRINT("q2 + q1 = ";)
    PRINTQ(T.)

    MULT_REAL(Q.,REAL->T.)
    PRINT("q * real =";)
    PRINTQ(T.)

! multiplication is not commutative
    MULT(Q1.,Q2.->R.)
    PRINT("q1 * q2=";)
    PRINTQ(R.)

    MULT(Q2.,Q1.->R.)
    PRINT("q2 * q1=";)
    PRINTQ(R.)
END PROGRAM

Euphoria

function norm(sequence q)
    return sqrt(power(q[1],2)+power(q[2],2)+power(q[3],2)+power(q[4],2))
end function

function conj(sequence q)
    q[2..4] = -q[2..4]
    return q
end function

function add(object q1, object q2)
    if atom(q1) != atom(q2) then
        if atom(q1) then
            q1 = {q1,0,0,0}
        else
            q2 = {q2,0,0,0}
        end if
    end if
    return q1+q2
end function

function mul(object q1, object q2)
    if sequence(q1) and sequence(q2) then
        return { q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3] - q1[4]*q2[4],
                 q1[1]*q2[2] + q1[2]*q2[1] + q1[3]*q2[4] - q1[4]*q2[3],
                 q1[1]*q2[3] - q1[2]*q2[4] + q1[3]*q2[1] + q1[4]*q2[2],
                 q1[1]*q2[4] + q1[2]*q2[3] - q1[3]*q2[2] + q1[4]*q2[1] }
    else
        return q1*q2
    end if
end function

function quats(sequence q)
    return sprintf("%g + %gi + %gj + %gk",q)
end function

constant
    q  = {1, 2, 3, 4},
    q1 = {2, 3, 4, 5},
    q2 = {5, 6, 7, 8},
    r  = 7

printf(1, "norm(q) = %g\n", norm(q))
printf(1, "-q = %s\n", {quats(-q)})
printf(1, "conj(q) = %s\n", {quats(conj(q))})
printf(1, "q + r = %s\n", {quats(add(q,r))})
printf(1, "q1 + q2 = %s\n", {quats(add(q1,q2))})
printf(1, "q1 * q2 = %s\n", {quats(mul(q1,q2))})
printf(1, "q2 * q1 = %s\n", {quats(mul(q2,q1))})
Output:
norm(q) = 5.47723
-q = -1 + -2i + -3j + -4k
conj(q) = 1 + -2i + -3j + -4k
q + r = 8 + 2i + 3j + 4k
q1 + q2 = 7 + 9i + 11j + 13k
q1 * q2 = -76 + 24i + 40j + 38k
q2 * q1 = -76 + 30i + 28j + 44k

F#

Mainly a

Translation of: C#

On the minus side we have no way to define a conversion to Quaternion from any suitable (numeric) type.

On the plus side we can avoid the stuff to make the equality structual (from the referential equality default) by just declaring it as an attribute to the type and let the compiler handle the details.

open System

[<Struct; StructuralEquality; NoComparison>]
type Quaternion(r : float, i : float, j : float, k : float) =
    member this.A = r
    member this.B = i
    member this.C = j
    member this.D = k

    new (f : float) = Quaternion(f, 0., 0., 0.)

    static member (~-) (q : Quaternion) = Quaternion(-q.A, -q.B, -q.C, -q.D)
 
    static member (+) (q1 : Quaternion, q2 : Quaternion) =
        Quaternion(q1.A + q2.A, q1.B + q2.B, q1.C + q2.C, q1.D + q2.D)
    static member (+) (q : Quaternion, r : float) = q + Quaternion(r)
    static member (+) (r : float, q: Quaternion) = Quaternion(r) + q
 
    static member (*) (q1 : Quaternion, q2 : Quaternion) =
        Quaternion(
            q1.A * q2.A - q1.B * q2.B - q1.C * q2.C - q1.D * q2.D,
            q1.A * q2.B + q1.B * q2.A + q1.C * q2.D - q1.D * q2.C,
            q1.A * q2.C - q1.B * q2.D + q1.C * q2.A + q1.D * q2.B,
            q1.A * q2.D + q1.B * q2.C - q1.C * q2.B + q1.D * q2.A)
    static member (*) (q : Quaternion, r : float) = q * Quaternion(r)
    static member (*) (r : float, q: Quaternion) = Quaternion(r) * q
 
    member this.Norm = Math.Sqrt(r * r + i * i + j * j + k * k)
 
    member this.Conjugate = Quaternion(r, -i, -j, -k)
 
    override this.ToString() = sprintf "Q(%f, %f, %f, %f)" r i j k

[<EntryPoint>]
let main argv =
    let q = Quaternion(1., 2., 3., 4.)
    let q1 = Quaternion(2., 3., 4., 5.)
    let q2 = Quaternion(3., 4., 5., 6.)
    let r = 7.
 
    printfn "q = %A" q
    printfn "q1 = %A" q1
    printfn "q2 = %A" q2
    printfn "r = %A" r
 
    printfn "q.Norm = %A" q.Norm
    printfn "q1.Norm = %A" q1.Norm
    printfn "q2.Norm = %A" q2.Norm
 
    printfn "-q = %A" -q
    printfn "q.Conjugate = %A" q.Conjugate
 
    printfn "q + r = %A" (q + (Quaternion r))
    printfn "q1 + q2 = %A" (q1 + q2)
    printfn "q2 + q1 = %A" (q2 + q1)
 
    printfn "q * r = %A" (q * r)
    printfn "q1 * q2 = %A" (q1 * q2)
    printfn "q2 * q1 = %A" (q2 * q1)
 
    printfn "q1*q2 %s q2*q1" (if (q1 * q2) = (q2 * q1) then "=" else "<>")
    printfn "q %s Q(1.,2.,3.,4.)" (if q = Quaternion(1., 2., 3., 4.) then "=" else "<>")
    0
Output:
q = Q(1.000000, 2.000000, 3.000000, 4.000000)
q1 = Q(2.000000, 3.000000, 4.000000, 5.000000)
q2 = Q(3.000000, 4.000000, 5.000000, 6.000000)
r = 7.0
q.Norm = 5.477225575
q1.Norm = 7.348469228
q2.Norm = 9.273618495
-q = Q(-1.000000, -2.000000, -3.000000, -4.000000)
q.Conjugate = Q(1.000000, -2.000000, -3.000000, -4.000000)
q + r = Q(8.000000, 2.000000, 3.000000, 4.000000)
q1 + q2 = Q(5.000000, 7.000000, 9.000000, 11.000000)
q2 + q1 = Q(5.000000, 7.000000, 9.000000, 11.000000)
q * r = Q(7.000000, 14.000000, 21.000000, 28.000000)
q1 * q2 = Q(-56.000000, 16.000000, 24.000000, 26.000000)
q2 * q1 = Q(-56.000000, 18.000000, 20.000000, 28.000000)
q1*q2 <> q2*q1
q = Q(1.,2.,3.,4.)

Factor

The math.quaternions vocabulary provides words for treating sequences like quaternions. norm and vneg come from the math.vectors vocabulary. Oddly, I wasn't able to find a word for adding a real to a quaternion, so I wrote one.

USING: generalizations io kernel locals math.quaternions
math.vectors prettyprint sequences ;
IN: rosetta-code.quaternion-type

: show ( quot -- )
    [ unparse 2 tail but-last "= " append write ] [ call . ] bi
    ; inline

: 2show ( quots -- )
    [ 2curry show ] map-compose [ call ] each ; inline

: q+n ( q n -- q+n ) n>q q+ ;

[let
    { 1 2 3 4 } 7 { 2 3 4 5 } { 3 4 5 6 } :> ( q r q1 q2 )
    q [ norm ]
    q [ vneg ]
    q [ qconjugate ]
    [ curry show ] 2tri@
    {
        [ q  r  [ q+n ] ]
        [ q  r  [ q*n ] ]
        [ q1 q2 [ q+  ] ]
        [ q1 q2 [ q*  ] ]
        [ q2 q1 [ q*  ] ]
    } 2show
]
Output:
{ 1 2 3 4 } norm = 5.477225575051661
{ 1 2 3 4 } vneg = { -1 -2 -3 -4 }
{ 1 2 3 4 } qconjugate = { 1 -2 -3 -4 }
{ 1 2 3 4 } 7 q+n = { 8 2 3 4 }
{ 1 2 3 4 } 7 q*n = { 7 14 21 28 }
{ 2 3 4 5 } { 3 4 5 6 } q+ = { 5 7 9 11 }
{ 2 3 4 5 } { 3 4 5 6 } q* = { -56 16 24 26 }
{ 3 4 5 6 } { 2 3 4 5 } q* = { -56 18 20 28 }

Forth

: quaternions  4 * floats ;

: qvariable create 1 quaternions allot ;

: q! ( a b c d q -- )
  dup 3 floats + f!  dup 2 floats + f!  dup float+ f!  f! ;

: qcopy ( src dest -- ) 1 quaternions move ;

: qnorm  ( q -- f )
  0e 4 0 do  dup f@ fdup f* f+  float+ loop drop fsqrt ;

: qf* ( q f -- )
  4 0 do dup f@ fover f* dup f!  float+ loop fdrop drop ;

: qnegate ( q -- )  -1e qf* ;

: qconj ( q -- )
  float+ 3 0 do dup f@ fnegate dup f!  float+ loop drop ;

: qf+ ( q f -- ) dup f@ f+ f! ;

: q+ ( q1 q2 -- )
  4 0 do over f@ dup f@ f+ dup f!  float+ swap float+ swap loop 2drop ;

\ access
: q.a             f@ ;
: q.b      float+ f@ ;
: q.c  2 floats + f@ ;
: q.d  3 floats + f@ ;

: q* ( dest q1 q2 -- )
  over q.a dup q.d f*  over q.b dup q.c f* f+  over q.c dup q.b f* f-  over q.d dup q.a f* f+
  over q.a dup q.c f*  over q.b dup q.d f* f-  over q.c dup q.a f* f+  over q.d dup q.b f* f+
  over q.a dup q.b f*  over q.b dup q.a f* f+  over q.c dup q.d f* f+  over q.d dup q.c f* f-
  over q.a dup q.a f*  over q.b dup q.b f* f-  over q.c dup q.c f* f-  over q.d dup q.d f* f-
  2drop  4 0 do dup f!  float+ loop  drop ;

: q= ( q1 q2 -- ? )
  4 0 do
    over f@ dup f@ f<> if 2drop false unloop exit then
    float+ swap float+
  loop
  2drop true ;

\ testing

: q. ( q -- )
  [char] ( emit space
  4 0 do dup f@ f.  float+ loop drop
  [char] ) emit space ;

qvariable q   1e 2e 3e 4e q  q!
qvariable q1  2e 3e 4e 5e q1 q!
create q2     3e f, 4e f, 5e f, 6e f,	\ by hand

qvariable tmp
qvariable m1
qvariable m2

q qnorm f.				\ 5.47722557505166
q tmp qcopy  tmp qnegate  tmp q.	\ ( -1. -2. -3. -4. )
q tmp qcopy  tmp qconj    tmp q.	\ ( 1. -2. -3. -4. )

q m1 qcopy  m1 7e qf+   m1 q.		\ ( 8. 2. 3. 4. )
q m2 qcopy  7e m2 qf+   m2 q.		\ ( 8. 2. 3. 4. )
m1 m2 q= .				\ -1  (true)

q2 tmp qcopy  q1 tmp q+   tmp q.	\ ( 5. 7. 9. 11. )

q m1 qcopy  m1 7e qf*     m1 q.		\ ( 7. 14. 21. 28. )
q m2 qcopy  7e m2 qf*     m2 q.		\ ( 7. 14. 21. 28. )
m1 m2 q= .				\ -1  (true)

m1 q1 q2 q*  m1 q.			\ ( -56. 16. 24. 26. )
m2 q2 q1 q*  m2 q.			\ ( -56. 18. 20. 28. )
m1 m2 q= .				\ 0  (false)

Fortran

Works with: Fortran version 90 and later
module Q_mod
  implicit none

  type quaternion
    real :: a, b, c, d
  end type

  public :: norm, neg, conj
  public :: operator (+)
  public :: operator (*)
  
  private ::  q_plus_q, q_plus_r, r_plus_q, &
              q_mult_q, q_mult_r, r_mult_q, &
              norm_q, neg_q, conj_q

  interface norm
    module procedure norm_q
  end interface

  interface neg
    module procedure neg_q
  end interface

  interface conj
    module procedure conj_q
  end interface

  interface operator (+)
    module procedure q_plus_q, q_plus_r, r_plus_q
  end interface

  interface operator (*)
    module procedure q_mult_q, q_mult_r, r_mult_q
  end interface

contains

function norm_q(x) result(res)
  real :: res
  type (quaternion), intent (in) :: x

  res = sqrt(x%a*x%a + x%b*x%b + x%c*x%c + x%d*x%d)
  
end function norm_q

function neg_q(x) result(res)
  type (quaternion) :: res
  type (quaternion), intent (in) :: x

  res%a = -x%a
  res%b = -x%b
  res%c = -x%c
  res%d = -x%d
  
end function neg_q

function conj_q(x) result(res)
  type (quaternion) :: res
  type (quaternion), intent (in) :: x

  res%a = x%a
  res%b = -x%b
  res%c = -x%c
  res%d = -x%d
  
end function conj_q

function q_plus_q(x, y) result (res)
  type (quaternion) :: res
  type (quaternion), intent (in) :: x, y
   
  res%a = x%a + y%a
  res%b = x%b + y%b
  res%c = x%c + y%c
  res%d = x%d + y%d
   
end function q_plus_q

function q_plus_r(x, r) result (res)
  type (quaternion) :: res 
  type (quaternion), intent (in) :: x
  real, intent(in) :: r
   
   res = x
   res%a = x%a + r
   
end function q_plus_r

function r_plus_q(r, x) result (res)
  type (quaternion) :: res 
  type (quaternion), intent (in) :: x
  real, intent(in) :: r
   
   res = x
   res%a = x%a + r
   
end function r_plus_q

function q_mult_q(x, y) result (res)
  type (quaternion) :: res 
  type (quaternion), intent (in) :: x, y
   
   res%a = x%a*y%a - x%b*y%b - x%c*y%c - x%d*y%d
   res%b = x%a*y%b + x%b*y%a + x%c*y%d - x%d*y%c
   res%c = x%a*y%c - x%b*y%d + x%c*y%a + x%d*y%b
   res%d = x%a*y%d + x%b*y%c - x%c*y%b + x%d*y%a
   
end function q_mult_q

function q_mult_r(x, r) result (res)
  type (quaternion) :: res 
  type (quaternion), intent (in) :: x
  real, intent(in) ::  r
   
   res%a = x%a*r
   res%b = x%b*r
   res%c = x%c*r  
   res%d = x%d*r 
   
end function q_mult_r

function r_mult_q(r, x) result (res)
  type (quaternion) :: res 
  type (quaternion), intent (in) :: x
  real, intent(in) ::  r
   
   res%a = x%a*r
   res%b = x%b*r
   res%c = x%c*r  
   res%d = x%d*r 
   
end function r_mult_q
end module Q_mod

program Quaternions
  use Q_mod
  implicit none

  real :: r = 7.0
  type(quaternion) :: q, q1, q2

  q  = quaternion(1, 2, 3, 4)
  q1 = quaternion(2, 3, 4, 5)
  q2 = quaternion(3, 4, 5, 6)

  write(*, "(a, 4f8.3)") "             q = ", q
  write(*, "(a, 4f8.3)") "            q1 = ", q1
  write(*, "(a, 4f8.3)") "            q2 = ", q2
  write(*, "(a, f8.3)")  "             r = ", r
  write(*, "(a, f8.3)")  "     Norm of q = ", norm(q) 
  write(*, "(a, 4f8.3)") " Negative of q = ", neg(q)
  write(*, "(a, 4f8.3)") "Conjugate of q = ", conj(q) 
  write(*, "(a, 4f8.3)") "         q + r = ", q + r
  write(*, "(a, 4f8.3)") "         r + q = ", r + q
  write(*, "(a, 4f8.3)") "       q1 + q2 = ", q1 + q2
  write(*, "(a, 4f8.3)") "         q * r = ", q * r  
  write(*, "(a, 4f8.3)") "         r * q = ", r * q 
  write(*, "(a, 4f8.3)") "       q1 * q2 = ", q1 * q2 
  write(*, "(a, 4f8.3)") "       q2 * q1 = ", q2 * q1 

end program
Output:
             q =    1.000   2.000   3.000   4.000
            q1 =    2.000   3.000   4.000   5.000
            q2 =    3.000   4.000   5.000   6.000
             r =    7.000
     Norm of q =    5.477
 Negative of q =   -1.000  -2.000  -3.000  -4.000
Conjugate of q =    1.000  -2.000  -3.000  -4.000
         q + r =    8.000   2.000   3.000   4.000
         r + q =    8.000   2.000   3.000   4.000
       q1 + q2 =    5.000   7.000   9.000  11.000
         q * r =    7.000  14.000  21.000  28.000
         r * q =    7.000  14.000  21.000  28.000
       q1 * q2 =  -56.000  16.000  24.000  26.000
       q2 * q1 =  -56.000  18.000  20.000  28.000

FreeBASIC

Dim Shared As Integer q(3)  = {1, 2, 3, 4}
Dim Shared As Integer q1(3) = {2, 3, 4, 5}
Dim Shared As Integer q2(3) = {3, 4, 5, 6}
Dim Shared As Integer i, r = 7, t(3)

Function q_norm(q() As Integer) As Double
    ' medida o valor absoluto de un cuaternión
    Dim As Double a = 0
    For i = 0 To 3 
        a += q(i)^2
    Next i
    Return Sqr(a)
End Function

Sub q_neg(q() As Integer) 
    For i = 0 To 3 
        q(i) *= -1
    Next i
End Sub

Sub q_conj(q() As Integer)
    ' conjugado de un cuaternión
    For i = 1 To 3 
        q(i) *= -1
    Next i
End Sub

Sub q_addreal(q() As Integer, r As Integer)
    q(0) += r 
End Sub

Sub q_add(q() As Integer, r() As Integer)
    ' adición entre cuaternios
    For i = 0 To 3 
        q(i) += r(i)
    Next i
End Sub

Sub q_mulreal(q() As Integer, r As Integer)
    For i = 0 To 3 
        q(i) *= r
    Next i
End Sub

Sub q_mul(q() As Integer, r() As Integer)
    ' producto entre cuaternios
    Dim As Integer m(3)
    m(0) = q(0)*r(0) - q(1)*r(1) - q(2)*r(2) - q(3)*r(3) 
    m(1) = q(0)*r(1) + q(1)*r(0) + q(2)*r(3) - q(3)*r(2) 
    m(2) = q(0)*r(2) - q(1)*r(3) + q(2)*r(0) + q(3)*r(1) 
    m(3) = q(0)*r(3) + q(1)*r(2) - q(2)*r(1) + q(3)*r(0) 
    For i = 0 To 3 : q(i) = m(i) : Next i
End Sub

Function q_show(q() As Integer) As String
    Dim As String a = "("
    For i = 0 To 3 
        a += Str(q(i)) + ", " 
    Next i
    Return Mid(a,1,Len(a)-2) + ")"
End Function

'--- Programa Principal ---
Print " q = "; q_show(q())
Print "q1 = "; q_show(q1())
Print "q2 = "; q_show(q2())
Print " r = "; r
Print "norm(q) ="; q_norm(q())
For i = 0 To 3 : t(i) = q(i) : Next i : q_neg(t())  : Print " neg(q) = "; q_show(t())
For i = 0 To 3 : t(i) = q(i) : Next i : q_conj(t()) : Print "conj(q) = "; q_show(t())
For i = 0 To 3 : t(i) = q(i) : Next i : q_addreal(t(),r) : Print " r + q  = "; q_show(t())
For i = 0 To 3 : t(i) = q1(i) : Next i : q_add(t(),q2()) : Print "q1 + q2 = "; q_show(t())
For i = 0 To 3 : t(i) = q2(i) : Next i : q_add(t(),q1()) : Print "q2 + q1 = "; q_show(t())
For i = 0 To 3 : t(i) = q(i) : Next i : q_mulreal(t(),r) : Print " r * q  = "; q_show(t())
For i = 0 To 3 : t(i) = q1(i) : Next i : q_mul(t(),q2()) : Print "q1 * q2 = "; q_show(t())
For i = 0 To 3 : t(i) = q2(i) : Next i : q_mul(t(),q1()) : Print "q2 * q1 = "; q_show(t())
End
Output:
 q = (1, 2, 3, 4)
q1 = (2, 3, 4, 5)
q2 = (3, 4, 5, 6)
r = 7
norm(q) = 5.477225575051661
neg(q)  = (-1, -2, -3, -4)
conj(q) = (1, -2, -3, -4)
 r + q  = (8, 2, 3, 4)
q1 + q2 = (5, 7, 9, 11)
q2 + q1 = (5, 7, 9, 11)
 r * q  = (7, 14, 21, 28)
q1 * q2 = (-56, 16, 24, 26)
q2 * q1 = (-56, 18, 20, 28)

GAP

# GAP has built-in support for quaternions

A := QuaternionAlgebra(Rationals);
# <algebra-with-one of dimension 4 over Rationals>

b := BasisVectors(Basis(A));
# [ e, i, j, k ]

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# Conjugate
ComplexConjugate(q);
# e+(-2)*i+(-3)*j+(-4)*k

# Division
1/q;
# (1/30)*e+(-1/15)*i+(-1/10)*j+(-2/15)*k

# Computing norm may be difficult, since the result would be in a quadratic field.
# Sqrt exists in GAP, but it is quite unusual: see ?E in GAP documentation, and the following example
Sqrt(5/3);
# 1/3*E(60)^7+1/3*E(60)^11-1/3*E(60)^19-1/3*E(60)^23-1/3*E(60)^31+1/3*E(60)^43-1/3*E(60)^47+1/3*E(60)^59

# However, the square of the norm is easy to compute
q*ComplexConjugate(q);
# (30)*e

q1 := [2, 3, 4, 5]*b;
# (2)*e+(3)*i+(4)*j+(5)*k

q2 := [3, 4, 5, 6]*b;
# (3)*e+(4)*i+(5)*j+(6)*k

q1*q2 - q2*q1;
# (-2)*i+(4)*j+(-2)*k

# Can't add directly to a rational, one must make a quaternion of it
r := 5/3*b[1];
# (5/3)*e
r + q;
# (8/3)*e+(2)*i+(3)*j+(4)*k

# For multiplication, no problem (we are in an algebra over rationals !)
r*q;
# (5/3)*e+(10/3)*i+(5)*j+(20/3)*k
5/3*q;
# (5/3)*e+(10/3)*i+(5)*j+(20/3)*k

# Negative
-q;
(-1)*e+(-2)*i+(-3)*j+(-4)*k


# While quaternions are built-in, you can define an algebra in GAP by specifying it's multiplication table.
# See tutorial, p. 60, and reference of the functions used below.

# A multiplication table of dimension 4.

T := EmptySCTable(4, 0);
SetEntrySCTable(T, 1, 1, [1, 1]);
SetEntrySCTable(T, 1, 2, [1, 2]);
SetEntrySCTable(T, 1, 3, [1, 3]);
SetEntrySCTable(T, 1, 4, [1, 4]);
SetEntrySCTable(T, 2, 1, [1, 2]);
SetEntrySCTable(T, 2, 2, [-1, 1]);
SetEntrySCTable(T, 2, 3, [1, 4]);
SetEntrySCTable(T, 2, 4, [-1, 3]);
SetEntrySCTable(T, 3, 1, [1, 3]);
SetEntrySCTable(T, 3, 2, [-1, 4]);
SetEntrySCTable(T, 3, 3, [-1, 1]);
SetEntrySCTable(T, 3, 4, [1, 2]);
SetEntrySCTable(T, 4, 1, [1, 4]);
SetEntrySCTable(T, 4, 2, [1, 3]);
SetEntrySCTable(T, 4, 3, [-1, 2]);
SetEntrySCTable(T, 4, 4, [-1, 1]);

A := AlgebraByStructureConstants(Rationals, T, ["e", "i", "j", "k"]);
b := GeneratorsOfAlgebra(A);

IsAssociative(A);
# true

IsCommutative(A);
# false

# Then, like above

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# However, as is, GAP does not know division or conjugate on this algebra.
# QuaternionAlgebra is useful as well for extensions of rationals,
# and this one _has_ conjugate and division, as seen previously.

# Try this on Q[z] where z is the square root of 5 (in GAP it's ER(5))
F := FieldByGenerators([ER(5)]);
A := QuaternionAlgebra(F);
b := GeneratorsOfAlgebra(A);

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# Conjugate and division

ComplexConjugate(q);
# e+(-2)*i+(-3)*j+(-4)*k

1/q;
# (1/30)*e+(-1/15)*i+(-1/10)*j+(-2/15)*k

Go

Conventions for method receiver, parameter, and return values modeled after Go's big number package. It provides flexibility without requiring unnecessary object creation. The test program creates only four quaternion objects, the three inputs and one more for an output. The three inputs are reused repeatedly without being modified. The output is also reused repeatedly, being overwritten for each operation.

package main

import (
    "fmt"
    "math"
)

type qtn struct {
    r, i, j, k float64
}

var (
    q  = &qtn{1, 2, 3, 4}
    q1 = &qtn{2, 3, 4, 5}
    q2 = &qtn{3, 4, 5, 6}

    r  float64 = 7
)

func main() {
    fmt.Println("Inputs")
    fmt.Println("q:", q)
    fmt.Println("q1:", q1)
    fmt.Println("q2:", q2)
    fmt.Println("r:", r)

    var qr qtn
    fmt.Println("\nFunctions")
    fmt.Println("q.norm():", q.norm())
    fmt.Println("neg(q):", qr.neg(q))
    fmt.Println("conj(q):", qr.conj(q))
    fmt.Println("addF(q, r):", qr.addF(q, r))
    fmt.Println("addQ(q1, q2):", qr.addQ(q1, q2))
    fmt.Println("mulF(q, r):", qr.mulF(q, r))
    fmt.Println("mulQ(q1, q2):", qr.mulQ(q1, q2))
    fmt.Println("mulQ(q2, q1):", qr.mulQ(q2, q1))
}

func (q *qtn) String() string {
    return fmt.Sprintf("(%g, %g, %g, %g)", q.r, q.i, q.j, q.k)
}

func (q *qtn) norm() float64 {
    return math.Sqrt(q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k)
}

func (z *qtn) neg(q *qtn) *qtn {
    z.r, z.i, z.j, z.k = -q.r, -q.i, -q.j, -q.k
    return z
}

func (z *qtn) conj(q *qtn) *qtn {
    z.r, z.i, z.j, z.k = q.r, -q.i, -q.j, -q.k
    return z
}

func (z *qtn) addF(q *qtn, r float64) *qtn {
    z.r, z.i, z.j, z.k = q.r+r, q.i, q.j, q.k
    return z
}

func (z *qtn) addQ(q1, q2 *qtn) *qtn {
    z.r, z.i, z.j, z.k = q1.r+q2.r, q1.i+q2.i, q1.j+q2.j, q1.k+q2.k
    return z
}

func (z *qtn) mulF(q *qtn, r float64) *qtn {
    z.r, z.i, z.j, z.k = q.r*r, q.i*r, q.j*r, q.k*r
    return z
}

func (z *qtn) mulQ(q1, q2 *qtn) *qtn {
    z.r, z.i, z.j, z.k =
        q1.r*q2.r-q1.i*q2.i-q1.j*q2.j-q1.k*q2.k,
        q1.r*q2.i+q1.i*q2.r+q1.j*q2.k-q1.k*q2.j,
        q1.r*q2.j-q1.i*q2.k+q1.j*q2.r+q1.k*q2.i,
        q1.r*q2.k+q1.i*q2.j-q1.j*q2.i+q1.k*q2.r
    return z
}
Output:
Inputs
q: (1, 2, 3, 4)
q1: (2, 3, 4, 5)
q2: (3, 4, 5, 6)
r: 7

Functions
q.norm(): 5.477225575051661
neg(q): (-1, -2, -3, -4)
conj(q): (1, -2, -3, -4)
addF(q, r): (8, 2, 3, 4)
addQ(q1, q2): (5, 7, 9, 11)
mulF(q, r): (7, 14, 21, 28)
mulQ(q1, q2): (-56, 16, 24, 26)
mulQ(q2, q1): (-56, 18, 20, 28)

Haskell

import Control.Monad (join)

data Quaternion a =
  Q a a a a
  deriving (Show, Eq)

realQ :: Quaternion a -> a
realQ (Q r _ _ _) = r

imagQ :: Quaternion a -> [a]
imagQ (Q _ i j k) = [i, j, k]

quaternionFromScalar :: (Num a) => a -> Quaternion a
quaternionFromScalar s = Q s 0 0 0

listFromQ :: Quaternion a -> [a]
listFromQ (Q a b c d) = [a, b, c, d]

quaternionFromList :: [a] -> Quaternion a
quaternionFromList [a, b, c, d] = Q a b c d

normQ :: (RealFloat a) => Quaternion a -> a
normQ = sqrt . sum . join (zipWith (*)) . listFromQ

conjQ :: (Num a) => Quaternion a -> Quaternion a
conjQ (Q a b c d) = Q a (-b) (-c) (-d)

instance (RealFloat a) => Num (Quaternion a) where
  (Q a b c d) + (Q p q r s) = Q (a + p) (b + q) (c + r) (d + s)
  (Q a b c d) - (Q p q r s) = Q (a - p) (b - q) (c - r) (d - s)
  (Q a b c d) * (Q p q r s) =
    Q
    (a * p - b * q - c * r - d * s)
    (a * q + b * p + c * s - d * r)
    (a * r - b * s + c * p + d * q)
    (a * s + b * r - c * q + d * p)
  negate (Q a b c d)        = Q (-a) (-b) (-c) (-d)
  abs q                     = quaternionFromScalar (normQ q)
  signum (Q 0 0 0 0)        = 0
  signum q@(Q a b c d)      = Q (a/n) (b/n) (c/n) (d/n) where n = normQ q
  fromInteger n             = quaternionFromScalar (fromInteger n)

main :: IO ()
main = do
  let q, q1, q2 :: Quaternion Double
      q  = Q 1 2 3 4
      q1 = Q 2 3 4 5
      q2 = Q 3 4 5 6
  print $ (Q 0 1 0 0) * (Q 0 0 1 0) * (Q 0 0 0 1) -- i*j*k; prints "Q (-1.0) 0.0 0.0 0.0"
  print $ q1 * q2                                 -- prints "Q (-56.0) 16.0 24.0 26.0"
  print $ q2 * q1                                 -- prints "Q (-56.0) 18.0 20.0 28.0"
  print $ q1 * q2 == q2 * q1                      -- prints "False"
  print $ imagQ q                                 -- prints "[2.0,3.0,4.0]"

Icon and Unicon

Using Unicon's class system.

class Quaternion(a, b, c, d)

  method norm ()
    return sqrt (a*a + b*b + c*c + d*d)
  end

  method negative ()
    return Quaternion(-a, -b, -c, -d)
  end

  method conjugate ()
    return Quaternion(a, -b, -c, -d)
  end

  method add (n) 
    if type(n) == "Quaternion__state" 
      then return Quaternion(a+n.a, b+n.b, c+n.c, d+n.d)
      else return Quaternion(a+n, b, c, d)
  end

  method multiply (n) 
    if type(n) == "Quaternion__state" 
      then return Quaternion(a*n.a - b*n.b - c*n.c - d*n.d,
                             a*n.b + b*n.a + c*n.d - d*n.c,
                             a*n.c - b*n.d + c*n.a + d*n.b,
                             a*n.d + b*n.c - c*n.b + d*n.a)
      else return Quaternion(a*n, b*n, c*n, d*n)
  end

  method sign (n)
    return if n >= 0 then "+" else "-"
  end

  method string () 
    return ("" || a || sign(b) || abs(b) || "i" || sign(c) || abs(c) || "j" || sign(d) || abs(d) || "k");
  end

  initially(a, b, c, d)
    self.a := if /a then 0 else a
    self.b := if /b then 0 else b
    self.c := if /c then 0 else c
    self.d := if /d then 0 else d
end

To test the above:

procedure main ()
  q := Quaternion (1,2,3,4)
  q1 := Quaternion (2,3,4,5)
  q2 := Quaternion (3,4,5,6)
  r := 7

  write ("The norm      of " || q.string() || " is " || q.norm ())
  write ("The negative  of " || q.string() || " is " || q.negative().string ())
  write ("The conjugate of " || q.string() || " is " || q.conjugate().string ())
  write ("Sum of " || q.string() || " and " || r || " is " || q.add(r).string ())
  write ("Sum of " || q.string() || " and " || q1.string() || " is " || q.add(q1).string ())
  write ("Product of " || q.string() || " and " || r || " is " || q.multiply(r).string ())
  write ("Product of " || q.string() || " and " || q1.string() || " is " || q.multiply(q1).string ())
  write ("q1*q2 = " || q1.multiply(q2).string ())
  write ("q2*q1 = " || q2.multiply(q1).string ())
end
Output:
The norm      of 1+2i+3j+4k is 5.477225575
The negative  of 1+2i+3j+4k is -1-2i-3j-4k
The conjugate of 1+2i+3j+4k is 1-2i-3j-4k
Sum of 1+2i+3j+4k and 7 is 8+2i+3j+4k
Sum of 1+2i+3j+4k and 2+3i+4j+5k is 3+5i+7j+9k
Product of 1+2i+3j+4k and 7 is 7+14i+21j+28k
Product of 1+2i+3j+4k and 2+3i+4j+5k is -36+6i+12j+12k
q1*q2 = -56+16i+24j+26k
q2*q1 = -56+18i+20j+28k

Idris

With dependent types we can implement the more general Cayley-Dickson construction. Here the dependent type CD n a is implemented. It depends on a natural number n, which is the number of iterations carried out, and the base type a. So the real numbers are just CD 0 Double, the complex numbers CD 1 Double and the quaternions CD 2 Double

module CayleyDickson

data CD : Nat -> Type -> Type where
    CDBase : a -> CD 0 a
    CDProd : CD n a -> CD n a -> CD (S n) a

pairTy : Nat -> Type -> Type
pairTy Z a = a
pairTy (S n) a = let b = pairTy n a in (b, b)

fromPair : (n : Nat) -> pairTy n a -> CD n a
fromPair Z x = CDBase x
fromPair (S m) (x, y) = CDProd (fromPair m x) $ fromPair m y

toPair : CD n a -> pairTy n a
toPair (CDBase x) = x
toPair (CDProd x v) = (toPair x, toPair v)

first : CD n a -> a
first (CDBase x) = x
first (CDProd x v) = first x
    
fromBase : Num a => (n : Nat) -> a -> CD n a
fromBase Z x = CDBase x
fromBase (S m) x = CDProd (fromBase m x) $ fromBase m 0

multSclr : Num a => CD n a -> a -> CD n a
multSclr (CDBase x) y = CDBase $ x * y
multSclr (CDProd x v) y = CDProd (multSclr x y) $ multSclr v y

divSclr : Fractional a => CD n a -> a -> CD n a
divSclr (CDBase x) y = CDBase $ x / y
divSclr (CDProd x v) y = CDProd (divSclr x y) $ divSclr v y

plusCD : Num a => CD n a -> CD n a -> CD n a
plusCD (CDBase x) (CDBase y) = CDBase $ x + y
plusCD (CDProd x v) (CDProd y w) = CDProd (plusCD x y) $ plusCD v w

negCD : Neg a => CD n a -> CD n a
negCD (CDBase x) = CDBase $ negate x
negCD (CDProd x v) = CDProd (negCD x) $ negCD v

minusCD : Neg a => CD n a -> CD n a -> CD n a
minusCD (CDBase x) (CDBase y) = CDBase $ x - y
minusCD (CDProd x v) (CDProd y w) = CDProd (minusCD x y) $ minusCD v w

conjCD : Neg a => CD n a -> CD n a
conjCD (CDBase x) = CDBase x
conjCD (CDProd x v) = CDProd (conjCD x) $ negCD v

multCD : Neg a => CD n a -> CD n a -> CD n a
multCD (CDBase x) (CDBase y) = CDBase $ x * y
multCD (CDProd x v) (CDProd y w) = CDProd (minusCD (multCD x y) (multCD (conjCD w) v)) $ plusCD (multCD w x) $ multCD v $ conjCD y

absSqrCD : Neg a => CD n a -> CD n a
absSqrCD x = multCD x $ conjCD x

sqrLnCD : Neg a => CD n a -> a
sqrLnCD = first . absSqrCD

recipCD : Neg a => Fractional a => CD n a -> CD n a
recipCD x = conjCD $ divSclr x $ sqrLnCD x

divCD : Neg a => Fractional a => CD n a -> CD n a -> CD n a
divCD x y = multCD x $ recipCD y

absCD : CD n Double -> Double
absCD x = sqrt $ sqrLnCD x

showComps : Show a => CD n a -> String
showComps (CDBase x) = show x
showComps (CDProd x v) = showComps x ++ ", " ++ showComps v

Eq a => Eq (CD n a) where
    (CDBase x) == (CDBase y) = x == y
    (CDProd x v) == (CDProd y w) = x == y && v == w

Show a => Show (CD n a) where
    show x = "(" ++ showComps x ++ ")"

Neg a => Num (CD n a) where
    (+) = plusCD
    (*) = multCD
    fromInteger m {n} = fromBase n $ fromInteger m

Neg a => Neg (CD n a) where
    negate = negCD
    (-) = minusCD

(Neg a, Fractional a) => Fractional (CD n a) where
    (/) = divCD
    recip = recipCD

Abs (CD n Double) where
    abs {n} = fromBase n . absCD

To test it:

import CayleyDickson

main : IO ()
main =
	do
		let q = fromPair 2 ((1, 2), (3, 4))
		let q1 = fromPair 2 ((2, 3), (4, 5))
		let q2 = fromPair 2 ((3, 4), (5, 6))
		printLn $ q1 * q2
		printLn $ q2 * q1
		printLn $ q1 * q2 == q2 * q1

J

Derived from the j wiki:

   NB. utilities
   ip=:   +/ .*             NB. inner product
   T=. (_1^#:0 10 9 12)*0 7 16 23 A.=i.4
   toQ=:  4&{."1 :[:        NB. real scalars -> quaternion

   NB. task
   norm=: %:@ip~@toQ        NB. | y
   neg=:  -&toQ             NB. - y  and  x - y
   conj=: 1 _1 _1 _1 * toQ  NB. + y
   add=:  +&toQ             NB. x + y
   mul=:  (ip T ip ])&toQ   NB. x * y

T is a rank 3 tensor which allows us to express quaternion product ab as the inner product ATB if A and B are 4 element vectors representing the quaternions a and b. (Note also that once we have defined mul we no longer need to retain the definition of T, so we define T using =. instead of =:). The value of T is probably more interesting than its definition, so:

   T
1  0  0  0
0  1  0  0
0  0  1  0
0  0  0  1

0 _1  0  0
1  0  0  0
0  0  0 _1
0  0  1  0

0  0 _1  0
0  0  0  1
1  0  0  0
0 _1  0  0

0  0  0 _1
0  0 _1  0
0  1  0  0
1  0  0  0

In other words, the last dimension of T corresponds to the structure of the right argument (columns, in the display of T), the first dimension of T corresponds to the structure of the left argument (tables, in the display of T) and the middle dimension of T corresponds to the structure of the result (rows, in the display of T).

Example use:

   q=: 1 2 3 4
   q1=: 2 3 4 5
   q2=: 3 4 5 6
   r=: 7
   
   norm q
5.47723
   neg q
_1 _2 _3 _4
   conj q
1 _2 _3 _4
   r add q
8 2 3 4
   q1 add q2
5 7 9 11
   r mul q
7 14 21 28
   q1 mul q2
_56 16 24 26
   q2 mul q1
_56 18 20 28

Finally, note that when quaternions are used to represent orientation or rotation, we are typically only interested in unit length quaternions. As this is the typical application for quaternions, you will sometimes see quaternion multiplication expressed using "simplifications" which are only valid for unit length quaternions. But note also that in many of those contexts you also need to normalize the quaternion length after multiplication.

(An exception to this need to normalize unit length quaternions after multiplication might be when quaternions are represented as an index into a geodesic grid. For example, a grid with 16x20 faces would have a total of 15 vertices for each face (5+4+3+2+1), 3 of those vertices would be from the original 20 vertices of the icosahedron, and 9 of those vertices (5+4+3-3) would be on the edge of the original face (and, thus, used for two faces), the remaining 3 vertices would be interior. This means we would have 170 vertices (20+(20*9%2)+20*3, which would allow a quaternion to be represented in a single byte index into a list of 170 quaternions, and would allow quaternion multiplication to be represented as a 29kbyte lookup table. In some contexts - where quaternion multiplication is needed in high volume for secondary or tertiary issues (where precision isn't vital), such low accuracy quaternions might be adequate or even an advantage...)

Java

public class Quaternion {
    private final double a, b, c, d;

    public Quaternion(double a, double b, double c, double d) {
        this.a = a;
        this.b = b;
        this.c = c;
        this.d = d;
    }
    public Quaternion(double r) {
        this(r, 0.0, 0.0, 0.0);
    }

    public double norm() {
        return Math.sqrt(a * a + b * b + c * c + d * d);
    }

    public Quaternion negative() {
        return new Quaternion(-a, -b, -c, -d);
    }

    public Quaternion conjugate() {
        return new Quaternion(a, -b, -c, -d);
    }

    public Quaternion add(double r) {
        return new Quaternion(a + r, b, c, d);
    }
    public static Quaternion add(Quaternion q, double r) {
        return q.add(r);
    }
    public static Quaternion add(double r, Quaternion q) {
        return q.add(r);
    }
    public Quaternion add(Quaternion q) {
        return new Quaternion(a + q.a, b + q.b, c + q.c, d + q.d);
    }
    public static Quaternion add(Quaternion q1, Quaternion q2) {
        return q1.add(q2);
    }

    public Quaternion times(double r) {
        return new Quaternion(a * r, b * r, c * r, d * r);
    }
    public static Quaternion times(Quaternion q, double r) {
        return q.times(r);
    }
    public static Quaternion times(double r, Quaternion q) {
        return q.times(r);
    }
    public Quaternion times(Quaternion q) {
        return new Quaternion(
            a * q.a - b * q.b - c * q.c - d * q.d,
            a * q.b + b * q.a + c * q.d - d * q.c,
            a * q.c - b * q.d + c * q.a + d * q.b,
            a * q.d + b * q.c - c * q.b + d * q.a
        );
    }
    public static Quaternion times(Quaternion q1, Quaternion q2) {
        return q1.times(q2);
    }

    @Override
    public boolean equals(Object obj) {
        if (!(obj instanceof Quaternion)) return false;
        final Quaternion other = (Quaternion) obj;
        if (Double.doubleToLongBits(this.a) != Double.doubleToLongBits(other.a)) return false;
        if (Double.doubleToLongBits(this.b) != Double.doubleToLongBits(other.b)) return false;
        if (Double.doubleToLongBits(this.c) != Double.doubleToLongBits(other.c)) return false;
        if (Double.doubleToLongBits(this.d) != Double.doubleToLongBits(other.d)) return false;
        return true;
    }
    @Override
    public String toString() {
        return String.format("%.2f + %.2fi + %.2fj + %.2fk", a, b, c, d).replaceAll("\\+ -", "- ");
    }

    public String toQuadruple() {
        return String.format("(%.2f, %.2f, %.2f, %.2f)", a, b, c, d);
    }

    public static void main(String[] args) {
        Quaternion q = new Quaternion(1.0, 2.0, 3.0, 4.0);
        Quaternion q1 = new Quaternion(2.0, 3.0, 4.0, 5.0);
        Quaternion q2 = new Quaternion(3.0, 4.0, 5.0, 6.0);
        double r = 7.0;
        System.out.format("q       = %s%n", q);
        System.out.format("q1      = %s%n", q1);
        System.out.format("q2      = %s%n", q2);
        System.out.format("r       = %.2f%n%n", r);
        System.out.format("\u2016q\u2016     = %.2f%n", q.norm());
        System.out.format("-q      = %s%n", q.negative());
        System.out.format("q*      = %s%n", q.conjugate());
        System.out.format("q + r   = %s%n", q.add(r));
        System.out.format("q1 + q2 = %s%n", q1.add(q2));
        System.out.format("q \u00d7 r   = %s%n", q.times(r));
        Quaternion q1q2 = q1.times(q2);
        Quaternion q2q1 = q2.times(q1);
        System.out.format("q1 \u00d7 q2 = %s%n", q1q2);
        System.out.format("q2 \u00d7 q1 = %s%n", q2q1);
        System.out.format("q1 \u00d7 q2 %s q2 \u00d7 q1%n", (q1q2.equals(q2q1) ? "=" : "\u2260"));
    }
}
Output:
q       = 1.00 + 2.00i + 3.00j + 4.00k
q1      = 2.00 + 3.00i + 4.00j + 5.00k
q2      = 3.00 + 4.00i + 5.00j + 6.00k
r       = 7.00

‖q‖     = 5.48
-q      = -1.00 - 2.00i - 3.00j - 4.00k
q*      = 1.00 - 2.00i - 3.00j - 4.00k
q + r   = 8.00 + 2.00i + 3.00j + 4.00k
q1 + q2 = 5.00 + 7.00i + 9.00j + 11.00k
q × r   = 7.00 + 14.00i + 21.00j + 28.00k
q1 × q2 = -56.00 + 16.00i + 24.00j + 26.00k
q2 × q1 = -56.00 + 18.00i + 20.00j + 28.00k
q1 × q2 ≠ q2 × q1

JavaScript

A basic implementation that runs on modern JS engines is given here. Quaternion.js is a more complete implementation having full operator support and various optimizations.

var Quaternion = (function() {
    // The Q() function takes an array argument and changes it
    // prototype so that it becomes a Quaternion instance.  This is
    // scoped only for prototype member access.
    function Q(a) {
	a.__proto__ = proto;
	return a;
    }

    // Actual constructor.  This constructor converts its arguments to
    // an array, then that array to a Quaternion instance, then
    // returns that instance.  (using "new" with this constructor is
    // optional)
    function Quaternion() {
	return Q(Array.prototype.slice.call(arguments, 0, 4));
    }

    // Prototype for all Quaternions
    const proto = {
	// Inherits from a 4-element Array
	__proto__ : [0,0,0,0],

	// Properties -- In addition to Array[0..3] access, we
	// also define matching a, b, c, and d properties
	get a() this[0],
	get b() this[1],
	get c() this[2],
	get d() this[3],

	// Methods
	norm : function() Math.sqrt(this.map(function(x) x*x).reduce(function(x,y) x+y)),
	negate : function() Q(this.map(function(x) -x)),
	conjugate : function() Q([ this[0] ].concat(this.slice(1).map(function(x) -x))),
	add : function(x) {
	    if ("number" === typeof x) {
		return Q([ this[0] + x ].concat(this.slice(1)));
	    } else {
		return Q(this.map(function(v,i) v+x[i]));
	    }
	},
	mul : function(r) {
	    var q = this;
	    if ("number" === typeof r) {
		return Q(q.map(function(e) e*r));
	    } else {
		return Q([ q[0] * r[0] - q[1] * r[1] - q[2] * r[2] - q[3] * r[3],
			   q[0] * r[1] + q[1] * r[0] + q[2] * r[3] - q[3] * r[2],
			   q[0] * r[2] - q[1] * r[3] + q[2] * r[0] + q[3] * r[1],
			   q[0] * r[3] + q[1] * r[2] - q[2] * r[1] + q[3] * r[0] ]);
	    }
	},
	equals : function(q) this.every(function(v,i) v === q[i]),
	toString : function() (this[0] + " + " + this[1] + "i + "+this[2] + "j + " + this[3] + "k").replace(/\+ -/g, '- ')
    };

    Quaternion.prototype = proto;
    return Quaternion;
})();

Task/Example Usage:

var q = Quaternion(1,2,3,4);
var q1 = Quaternion(2,3,4,5);
var q2 = Quaternion(3,4,5,6);
var r = 7;

console.log("q = "+q);
console.log("q1 = "+q1);
console.log("q2 = "+q2);
console.log("r = "+r);
console.log("1. q.norm() = "+q.norm());
console.log("2. q.negate() = "+q.negate());
console.log("3. q.conjugate() = "+q.conjugate());
console.log("4. q.add(r) = "+q.add(r));
console.log("5. q1.add(q2) = "+q1.add(q2));
console.log("6. q.mul(r) = "+q.mul(r));
console.log("7.a. q1.mul(q2) = "+q1.mul(q2));
console.log("7.b. q2.mul(q1) = "+q2.mul(q1));
console.log("8. q1.mul(q2) " + (q1.mul(q2).equals(q2.mul(q1)) ? "==" : "!=") + " q2.mul(q1)");
Output:
q = 1 + 2i + 3j + 4k
q1 = 2 + 3i + 4j + 5k
q2 = 3 + 4i + 5j + 6k
r = 7
1. q.norm() = 5.477225575051661
2. q.negate() = -1 - 2i - 3j - 4k
3. q.conjugate() = 1 - 2i - 3j - 4k
4. q.add(r) = 8 + 2i + 3j + 4k
5. q1.add(q2) = 5 + 7i + 9j + 11k
6. q.mul(r) = 7 + 14i + 21j + 28k
7.a. q1.mul(q2) = -56 + 16i + 24j + 26k
7.b. q2.mul(q1) = -56 + 18i + 20j + 28k
8. q1.mul(q2) != q2.mul(q1)

jq

Program file: quaternion.jq

def Quaternion(q0;q1;q2;q3): { "q0": q0, "q1": q1, "q2": q2, "q3": q3, "type": "Quaternion" };

# promotion of a real number to a quaternion 
def Quaternion(r): if (r|type) == "number" then Quaternion(r;0;0;0) else r end;

# thoroughly recursive pretty-print
def pp:

  def signage: if . >= 0 then "+ \(.)" else  "- \(-.)" end;

  if type == "object" then
     if .type == "Quaternion" then
       "\(.q0) \(.q1|signage)i \(.q2|signage)j \(.q3|signage)k"
     else with_entries( {key, "value" : (.value|pp)} )
     end
  elif type == "array" then map(pp)
  else . 
  end ;

def real(z): Quaternion(z).q0;

# Note: imag(z) returns the "i" component only,
# reflecting the embedding of the complex numbers within the quaternions:
def imag(z): Quaternion(z).q1;

def conj(z): Quaternion(z) | Quaternion(.q0; -(.q1); -(.q2); -(.q3));

def abs2(z): Quaternion(z) | .q0 * .q0 + .q1*.q1 + .q2*.q2 + .q3*.q3;

def abs(z): abs2(z) | sqrt;

def negate(z): Quaternion(z) | Quaternion(-.q0; -.q1; -.q2; -.q3);

# z + w
def plus(z; w):
  def plusq(z;w): Quaternion(z.q0 + w.q0; z.q1 + w.q1;
                             z.q2 + w.q2; z.q3 + w.q3);
  plusq( Quaternion(z); Quaternion(w) );

# z - w
def minus(z; w):
  def minusq(z;w): Quaternion(z.q0 - w.q0; z.q1 - w.q1;
                              z.q2 - w.q2; z.q3 - w.q3);
  minusq( Quaternion(z); Quaternion(w) );

# *
def times(z; w):
  def timesq(z; w):
       Quaternion(z.q0*w.q0 - z.q1*w.q1 - z.q2*w.q2 - z.q3*w.q3;
                  z.q0*w.q1 + z.q1*w.q0 + z.q2*w.q3 - z.q3*w.q2;
                  z.q0*w.q2 - z.q1*w.q3 + z.q2*w.q0 + z.q3*w.q1;
                  z.q0*w.q3 + z.q1*w.q2 - z.q2*w.q1 + z.q3*w.q0);
  timesq( Quaternion(z); Quaternion(w) );

# (z/w)
def div(z; w):
  if (w|type) == "number" then Quaternion(z.q0/w; z.q1/w; z.q2/w; z.q3/w)
  else times(z; inv(w))
  end;

def inv(z): div(conj(z); abs2(z));


# Example usage and output:

def say(msg; e): "\(msg) => \(e|pp)";

def demo:
  say( "Quaternion(1;0;0;0)"; Quaternion(1;0;0;0)),
  (Quaternion (1; 2; 3; 4) as $q 
  | Quaternion(2; 3; 4; 5) as $q1
  | Quaternion(3; 4; 5; 6) as $q2
  | 7 as $r
  | say( "abs($q)";        abs($q) ),   # norm
    say( "negate($q)";     negate($q) ),
    say( "conj($q)";       conj($q) ),
    "",
    say( "plus($r; $q)";   plus($r; $q)),
    say( "plus($q; $r)";   plus($q; $r)),
    "",
    say( "plus($q1; $q2 )"; plus($q1; $q2)),
    "",
    say( "times($r;$q)";    times($r;$q)),
    say( "times($q;$r)";    times($q;$r)),
    "",
    say( "times($q1;$q2)";  times($q1;$q2)),
    say( "times($q2; $q1)"; times($q2; $q1)),
    say( "times($q1; $q2) != times($q2; $q1)";
         times($q1; $q2) != times($q2; $q1) )
    ) ;

demo

Example usage and output:

# jq -c -n -R -f quaternion.jq
Quaternion(1;0;0;0) => 1 + 0i + 0j + 0k
abs($q) => 5.477225575051661
negate($q) => -1 - 2i - 3j + -4k
conj($q) => 1 - 2i - 3j - 4k

plus($r; $q) => 8 + 2i + 3j + 4k
plus($q; $r) => 8 + 2i + 3j + 4k

plus($q1; $q2 ) => 5 + 7i + 9j + 11k

times($r;$q) => 7 + 14i + 21j + 28k
times($q;$r) => 7 + 14i + 21j + 28k

times($q1;$q2) => -56 + 16i + 24j + 26k
times($q2; $q1) => -56 + 18i + 20j + 28k
times($q1; $q2) != times($q2; $q1) => true

Julia

Quaternion.jl has a more complete implementation. This is derived from the quaternion example file included with Julia 0.2, which implements a quaternion type complete with arithmetic, type conversions / promotion rules, polymorphism over arbitrary real numeric types, and pretty-printing.

import Base: convert, promote_rule, show, conj, abs, +, -, *

immutable Quaternion{T<:Real} <: Number
    q0::T
    q1::T
    q2::T
    q3::T
end

Quaternion(q0::Real,q1::Real,q2::Real,q3::Real) = Quaternion(promote(q0,q1,q2,q3)...)

convert{T}(::Type{Quaternion{T}}, x::Real) =
    Quaternion(convert(T,x), zero(T), zero(T), zero(T))
convert{T}(::Type{Quaternion{T}}, z::Complex) =
    Quaternion(convert(T,real(z)), convert(T,imag(z)), zero(T), zero(T))
convert{T}(::Type{Quaternion{T}}, z::Quaternion) =
    Quaternion(convert(T,z.q0), convert(T,z.q1), convert(T,z.q2), convert(T,z.q3))

promote_rule{T,S}(::Type{Complex{T}}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}
promote_rule{T<:Real,S}(::Type{T}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}
promote_rule{T,S}(::Type{Quaternion{T}}, ::Type{Quaternion{S}}) = Quaternion{promote_type(T,S)}

function show(io::IO, z::Quaternion)
    pm(x) = x <	0 ? " - $(-x)" : " + $x"
    print(io, z.q0, pm(z.q1), "i", pm(z.q2), "j", pm(z.q3), "k")
end

conj(z::Quaternion) = Quaternion(z.q0, -z.q1, -z.q2, -z.q3)
abs(z::Quaternion) = sqrt(z.q0*z.q0 + z.q1*z.q1 + z.q2*z.q2 + z.q3*z.q3)

(-)(z::Quaternion) = Quaternion(-z.q0, -z.q1, -z.q2, -z.q3)

(+)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 + w.q0, z.q1 + w.q1,
                                               z.q2 + w.q2, z.q3 + w.q3)
(-)(z::Quaternion, w::Quaternion) = Quaternion(z.q0 - w.q0, z.q1 - w.q1,
                                               z.q2 - w.q2, z.q3 - w.q3)
(*)(z::Quaternion, w::Quaternion) = Quaternion(z.q0*w.q0 - z.q1*w.q1 - z.q2*w.q2 - z.q3*w.q3,
                                               z.q0*w.q1 + z.q1*w.q0 + z.q2*w.q3 - z.q3*w.q2,
                                               z.q0*w.q2 - z.q1*w.q3 + z.q2*w.q0 + z.q3*w.q1,
                                               z.q0*w.q3 + z.q1*w.q2 - z.q2*w.q1 + z.q3*w.q0)

Example usage and output:

julia> q = Quaternion(1,0,0,0)
julia> q  = Quaternion (1, 2, 3, 4)
       q1 = Quaternion(2, 3, 4, 5)
       q2 = Quaternion(3, 4, 5, 6)
       r = 7.

julia> norm(q)
5.477225575051661

julia> -q
-1 - 2i - 3j - 4k

julia> conj(q)
1 - 2i - 3j - 4k

julia> r + q, q + r
(8.0 + 2.0i + 3.0j + 4.0k,8.0 + 2.0i + 3.0j + 4.0k)

julia> q1 + q2
5 + 7i + 9j + 11k

julia> r*q, q*r
(7.0 + 14.0i + 21.0j + 28.0k,7.0 + 14.0i + 21.0j + 28.0k)

julia> q1*q2, q2*q1, q1*q2 != q2*q1
(-56 + 16i + 24j + 26k,-56 + 18i + 20j + 28k,true)

Kotlin

// version 1.1.2

data class Quaternion(val a: Double, val b: Double, val c: Double, val d: Double) {
    operator fun plus(other: Quaternion): Quaternion {
        return Quaternion (this.a + other.a, this.b + other.b,
                           this.c + other.c, this.d + other.d)
    }

    operator fun plus(r: Double) = Quaternion(a + r, b, c, d)

    operator fun times(other: Quaternion): Quaternion {
        return Quaternion(
            this.a * other.a - this.b * other.b - this.c * other.c - this.d * other.d,
            this.a * other.b + this.b * other.a + this.c * other.d - this.d * other.c,
            this.a * other.c - this.b * other.d + this.c * other.a + this.d * other.b,
            this.a * other.d + this.b * other.c - this.c * other.b + this.d * other.a
        )
    }

    operator fun times(r: Double) = Quaternion(a * r, b * r, c * r, d * r)

    operator fun unaryMinus() =  Quaternion(-a, -b, -c, -d)

    fun conj() = Quaternion(a, -b, -c, -d)

    fun norm() = Math.sqrt(a * a + b * b + c * c + d * d)

    override fun toString() = "($a, $b, $c, $d)"
}

// extension functions for Double type
operator fun Double.plus(q: Quaternion) = q + this
operator fun Double.times(q: Quaternion) = q * this

fun main(args: Array<String>) {
    val q  = Quaternion(1.0, 2.0, 3.0, 4.0)
    val q1 = Quaternion(2.0, 3.0, 4.0, 5.0)
    val q2 = Quaternion(3.0, 4.0, 5.0, 6.0)
    val r  = 7.0
    println("q  = $q")
    println("q1 = $q1")
    println("q2 = $q2")
    println("r  = $r\n")
    println("norm(q) = ${"%f".format(q.norm())}")
    println("-q      = ${-q}")
    println("conj(q) = ${q.conj()}\n")
    println("r  + q  = ${r + q}")
    println("q  + r  = ${q + r}")
    println("q1 + q2 = ${q1 + q2}\n")
    println("r  * q  = ${r * q}")
    println("q  * r  = ${q * r}")
    val q3 = q1 * q2
    val q4 = q2 * q1
    println("q1 * q2 = $q3")
    println("q2 * q1 = $q4\n")
    println("q1 * q2 != q2 * q1 = ${q3 != q4}")
}
Output:
q  = (1.0, 2.0, 3.0, 4.0)
q1 = (2.0, 3.0, 4.0, 5.0)
q2 = (3.0, 4.0, 5.0, 6.0)
r  = 7.0

norm(q) = 5.477226
-q      = (-1.0, -2.0, -3.0, -4.0)
conj(q) = (1.0, -2.0, -3.0, -4.0)

r  + q  = (8.0, 2.0, 3.0, 4.0)
q  + r  = (8.0, 2.0, 3.0, 4.0)
q1 + q2 = (5.0, 7.0, 9.0, 11.0)

r  * q  = (7.0, 14.0, 21.0, 28.0)
q  * r  = (7.0, 14.0, 21.0, 28.0)
q1 * q2 = (-56.0, 16.0, 24.0, 26.0)
q2 * q1 = (-56.0, 18.0, 20.0, 28.0)

q1 * q2 != q2 * q1 = true

Liberty BASIC

Quaternions saved as a space-separated string of four numbers.

 q$ = q$( 1 , 2 , 3 , 4 )
q1$ = q$( 2 , 3 , 4 , 5 )
q2$ = q$( 3 , 4 , 5 , 6 )

real = 7

print "q = "  ;  q$
print "q1 = " ; q1$
print "q2 = " ; q2$ 

print "real = " ; real

print "length /norm q  = " ; length( q$ )               '   =norm                        norm of q
print "negative (-q1)  = " ; negative$( q1$ )           '   =negative                    negated q1
print "conjugate q     = " ; conjugate$( q$ )           '   conjugate                    conjugate q
print "real + q        = " ; add1$( q$ , real )         '   real +quaternion             real +q
print "q + q2          = " ; add2$( q$ , q2$ )          '   sum two quaternions          q +q2
print "real * q        = " ; multiply1$( q$ , real )    '   real *quaternion             real *q
print "q1 * q2         = " ; multiply2$( q1$ , q2$ )    '   product of two quaternions   q1 & q2
print "q2 * q1         = " ; multiply2$( q2$ , q1$ )    '   show q1 *q2 <> q2 *q1

end

function q$( r , i , j , k )
  q$ = str$( r); " "; str$( i); " "; str$( j); " "; str$( k)
end function

function length( q$ )
  r = val( word$( q$ , 1 ) )
  i = val( word$( q$ , 2 ) )
  j = val( word$( q$ , 3 ) )
  k = val( word$( q$ , 4 ) )
  length =sqr( r^2 +i^2 +j^2 +k^2)
end function

function multiply1$( q$ , d )
  r = val( word$( q$ , 1 ) )
  i = val( word$( q$ , 2 ) )
  j = val( word$( q$ , 3 ) )
  k = val( word$( q$ , 4 ) )
  multiply1$ =q$( r*d, i*d, j*d, k*d)
end function

function multiply2$( q$ , b$ )
  ar = val( word$( q$ , 1 ) )   'a1
  ai = val( word$( q$ , 2 ) )   'b1
  aj = val( word$( q$ , 3 ) )   'c1
  ak = val( word$( q$ , 4 ) )   'd1

  br = val( word$( b$ , 1 ) )   'a2
  bi = val( word$( b$ , 2 ) )   'b2
  bj = val( word$( b$ , 3 ) )   'c2
  bk = val( word$( b$ , 4 ) )   'd2

  multiply2$ =q$( _
  ar *br_
  +( 0 -ai) *bi_
  +( 0 -aj) *bj_
  +( 0 -ak) *bk _
  ,_
  ar *bi_
  +ai *br_
  +aj *bk_
  +( 0 -ak) *bj_
  ,_
  ar *bj_
  +( 0 -ai) *bk_
  +aj *br_
  +ak *bi_
  ,_
  ar *bk_
  +ai *bj_
  +( 0 -aj) *bi_
  +ak *br )
end function

function negative$( q$ )
  r = val( word$( q$ , 1 ) )
  i = val( word$( q$ , 2 ) )
  j = val( word$( q$ , 3 ) )
  k = val( word$( q$ , 4 ) )
  negative$ =q$( 0-r, 0-i, 0-j, 0-k)
end function

function conjugate$( q$ )
  r = val( word$( q$ , 1 ) )
  i = val( word$( q$ , 2 ) )
  j = val( word$( q$ , 3 ) )
  k = val( word$( q$ , 4 ) )
  conjugate$ =q$( r, 0-i, 0-j, 0-k)
end function

function add1$( q$ , real )
  r = val( word$( q$ , 1 ) )
  i = val( word$( q$ , 2 ) )
  j = val( word$( q$ , 3 ) )
  k = val( word$( q$ , 4 ) )
  add1$ =q$( r +real, i, j, k)
end function

function add2$( q$ , b$ )
  ar = val( word$( q$ , 1 ) )
  ai = val( word$( q$ , 2 ) )
  aj = val( word$( q$ , 3 ) )
  ak = val( word$( q$ , 4 ) )
  br = val( word$( b$ , 1 ) )
  bi = val( word$( b$ , 2 ) )
  bj = val( word$( b$ , 3 ) )
  bk = val( word$( b$ , 4 ) )
  add2$ =q$( ar +br, ai +bi, aj +bj, ak +bk)
end function

Lua

Quaternion = {}

function Quaternion.new( a, b, c, d )
    local q = { a = a or 1, b = b or 0, c = c or 0, d = d or 0 }

    local metatab = {}
    setmetatable( q, metatab )
    metatab.__add = Quaternion.add
    metatab.__sub = Quaternion.sub
    metatab.__unm = Quaternion.unm
    metatab.__mul = Quaternion.mul

    return q
end

function Quaternion.add( p, q )
    if type( p ) == "number" then
	return Quaternion.new( p+q.a, q.b, q.c, q.d )
    elseif type( q ) == "number" then
	return Quaternion.new( p.a+q, p.b, p.c, p.d )
    else
	return Quaternion.new( p.a+q.a, p.b+q.b, p.c+q.c, p.d+q.d )
    end
end

function Quaternion.sub( p, q )
    if type( p ) == "number" then
	return Quaternion.new( p-q.a, q.b, q.c, q.d )
    elseif type( q ) == "number" then
	return Quaternion.new( p.a-q, p.b, p.c, p.d )
    else
	return Quaternion.new( p.a-q.a, p.b-q.b, p.c-q.c, p.d-q.d )
    end
end

function Quaternion.unm( p )
    return Quaternion.new( -p.a, -p.b, -p.c, -p.d )
end

function Quaternion.mul( p, q )
    if type( p ) == "number" then
	return Quaternion.new( p*q.a, p*q.b, p*q.c, p*q.d )
    elseif type( q ) == "number" then
	return Quaternion.new( p.a*q, p.b*q, p.c*q, p.d*q )
    else
	return Quaternion.new( p.a*q.a - p.b*q.b - p.c*q.c - p.d*q.d,
                               p.a*q.b + p.b*q.a + p.c*q.d - p.d*q.c,
 			       p.a*q.c - p.b*q.d + p.c*q.a + p.d*q.b,
			       p.a*q.d + p.b*q.c - p.c*q.b + p.d*q.a )
    end
end

function Quaternion.conj( p )
    return Quaternion.new( p.a, -p.b, -p.c, -p.d )
end

function Quaternion.norm( p )
    return math.sqrt( p.a^2 + p.b^2 + p.c^2 + p.d^2 )
end

function Quaternion.print( p )
    print( string.format( "%f + %fi + %fj + %fk\n", p.a, p.b, p.c, p.d ) )
end

Examples:

q1 = Quaternion.new( 1, 2, 3, 4 )
q2 = Quaternion.new( 5, 6, 7, 8 )
r  = 12

print( "norm(q1) = ", Quaternion.norm( q1 ) )
io.write( "-q1 = " ); Quaternion.print( -q1 )
io.write( "conj(q1) = " ); Quaternion.print( Quaternion.conj( q1 ) )
io.write( "r+q1 = " ); Quaternion.print( r+q1 )
io.write( "q1+r = " ); Quaternion.print( q1+r )
io.write( "r*q1 = " ); Quaternion.print( r*q1 )
io.write( "q1*r = " ); Quaternion.print( q1*r )
io.write( "q1*q2 = " ); Quaternion.print( q1*q2 )
io.write( "q2*q1 = " ); Quaternion.print( q2*q1 )
Output:
norm(q1) = 5.4772255750517
-q1 = -1.000000 -2.000000i -3.000000j -4.000000k
conj(q1) = 1.000000 -2.000000i -3.000000j -4.000000k
r+q1 = 13.000000 + 2.000000i + 3.000000j + 4.000000k
q1+r = 13.000000 + 2.000000i + 3.000000j + 4.000000k
r*q1 = 12.000000 + 24.000000i + 36.000000j + 48.000000k
q1*r = 12.000000 + 24.000000i + 36.000000j + 48.000000k
q1*q2 = -60.000000 + 12.000000i + 30.000000j + 24.000000k
q2*q1 = -60.000000 + 20.000000i + 14.000000j + 32.000000k

M2000 Interpreter

We can define Quaternions using a class, using operators for specific tasks, as negate, add, multiplication and equality with rounding to 13 decimal place (thats what doing "==" operator for doubles)

Module CheckIt {
      class Quaternion {
            \\ by default are double
            a,b,c,d
            Property ToString$ {
                  Value {
                      link parent a,b,c, d to a,b,c,d  
                       value$=format$("{0} + {1}i + {2}j + {3}k",a,b,c,d)
                  }
            }
            Property Norm { Value}
            Operator "==" {
                  read n
                  push .a==n.a and .b==n.b and .c==n.c and .d==n.d
            }
            Module CalcNorm {
                  .[Norm]<=sqrt(.a**2+.b**2+.c**2+.d**2)
             }
            Operator Unary {
                  .a-! : .b-! : .c-! :.d-!
            }
            Function Conj {
                  q=this 
                  for q {
                         .b-! : .c-! :.d-!
                  }
                  =q
            }
            Function Add {
                  q=this
                  for q {
                         .a+=Number : .CalcNorm
                  }
                  =q
            }
            Operator "+"  {
                  Read q2
                  For this, q2 {
                        .a+=..a :.b+=..b:.c+=..c:.d+=..d
                        .CalcNorm
                  }
            }
            Function Mul(r)  {
                  q=this 
                  for q {
                        .a*=r:.b*=r:.c*=r:.d*=r:.CalcNorm
                  }
                  =q
            }
            Operator "*" {
                  Read q2
                  For This, q2 {
                        Push .a*..a-.b*..b-.c*..c-.d*..d
                        Push .a*..b+.b*..a+.c*..d-.d*..c
                        Push .a*..c-.b*..d+.c*..a+.d*..b
                        .d<=.a*..d+.b*..c-.c*..b+.d*..a
                        Read .c, .b, .a
                        .CalcNorm 
                  }
            }      
            class:
            module Quaternion {
                  if match("NNNN") then {
                        Read .a,.b,.c,.d
                       .CalcNorm
                  }
            }
      }
      \\ variables
      r=7
      q=Quaternion(1,2,3,4)
      q1=Quaternion(2,3,4,5)
      q2=Quaternion(3,4,5,6)
      
      \\ perform negate, conjugate, multiply by real, add a real, multiply quanterions, multiply in reverse order
      qneg=-q
      qconj=q.conj()
      qmul=q.Mul(r)
      qadd=q.Add(r)
      q1q2=q1*q2
      q2q1=q2*q1
      
      Print "q = ";q.ToString$
      Print "Normal q = ";q.Norm
      Print "Neg q = ";qneg.ToString$
      Print "Conj q = ";qconj.ToString$
      Print "Mul q 7 = ";qmul.ToString$
      Print "Add q 7 = ";qadd.ToString$
      Print "q1 = ";q1.ToString$
      Print "q2 = ";q2.ToString$
      Print "q1 * q2 = ";q1q2.ToString$
      Print "q2 * q1 = ";q2q1.ToString$
      Print q1==q1   ' true
      Print q1q2==q2q1 ' false
      \\ multiplication and equality in one expression
      Print q1 * q2 == q2 * q1  ' false
      Print q1 * q2 == q1 * q2  ' true
}
CheckIt
Output:
q = 1 + 2i + 3j + 4k
Normal q = 5.47722557505166
Neg q = -1 + -2i + -3j + -4k
Conj q = 1 + -2i + -3j + -4k            
Mul q 7 = 7 + 14i + 21j + 28k
Add q 7 = 8 + 2i + 3j + 4k
q1 = 2 + 3i + 4j + 5k
q2 = 3 + 4i + 5j + 6k
q1 * q2 = -56 + 16i + 24j + 26k
q2 * q1 = -56 + 18i + 20j + 28k
     True
    False
    false
     True

Maple

with(ArrayTools);

module Quaternion()
 option object;
 local real := 0;
 local i := 0;
 local j := 0;
 local k := 0;

 export getReal::static := proc(self::Quaternion, $)
  return self:-real;
 end proc;

 export getI::static := proc(self::Quaternion, $)
  return self:-i;
 end proc;

 export getJ::static := proc(self::Quaternion, $)
  return self:-j;
 end proc;

 export getK::static := proc(self::Quaternion, $)
  return self:-k;
 end proc;

 export Norm::static := proc(self::Quaternion, $)
  return sqrt(self:-real^2 + self:-i^2 + self:-j^2 + self:-k^2);
 end proc;

 # NegativeQuaternion returns the additive inverse of the quaternion
 export NegativeQuaternion::static := proc(self::Quaternion, $)
  return Quaternion(- self:-real, - self:-i, - self:-j, - self:-k);
 end proc;

 export Conjugate::static := proc(self::Quaternion, $)
  return Quaternion(self:-real, - self:-i, - self:-j, - self:-k);
 end proc;

 # quaternion addition
 export `+`::static := overload ([
  proc(self::Quaternion, x::Quaternion) option overload;
   return Quaternion(self:-real + getReal(x), self:-i + getI(x), self:-j + getJ(x), self:-k + getK(x));
  end proc,
  proc(self::Quaternion, x::algebraic) option overload;
   return Quaternion(self:-real + x, self:-i, self:-j, self:-k);
  end proc,
  proc(x::algebraic, self::Quaternion) option overload;
   return Quaternion(x + self:-real, self:-i, self:-j, self:-k);
  end
 ]);

 # convert quaternion to additive inverse
 export `-`::static := overload([
  proc(self::Quaternion) option overload;
   return Quaternion(-self:-real, -self:-i, -self:-j, -self:-k);
  end
 ]);

 # quaternion multiplication is non-abelian so the `.` operator needs to be used
 export `.`::static := overload([
  proc(self::Quaternion, x::Quaternion) option overload;
  return Quaternion(self:-real * getReal(x) - self:-i * getI(x) - self:-j * getJ(x) - self:-k * getK(x),
                    self:-real * getI(x) + self:-i * getReal(x) + self:-j * getK(x) - self:-k * getJ(x),
                    self:-real * getJ(x) + self:-j * getReal(x) - self:-i * getK(x) + self:-k * getI(x),
                    self:-real * getK(x) + self:-k * getReal(x) + self:-i * getJ(x) - self:-j * getI(x));
  end proc,
  proc(self::Quaternion, x::algebraic) option overload;
   return Quaternion(self:-real * x, self:-i * x, self:-j * x, self:-k * x);
  end proc,
  proc(x::algebraic, self::Quaternion) option overload;
   return Quaternion(self:-real * x, self:-i * x, self:-j * x, self:-k * x);
  end
 ]);

 # redirect division to `.` operator
 export `*`::static := overload([
  proc(self::Quaternion, x::Quaternion) option overload;
   use `*` = `.` in return self * x; end use
  end proc,
  proc(self::Quaternion, x::algebraic) option overload;
   use `*` = `.` in return x * self; end use
  end proc,
  proc(x::algebraic, self::Quaternion) option overload;
   use `*` = `.` in return x * self; end use
  end
 ]);

 # convert quaternion to multiplicative inverse
 export `/`::static := overload([
  proc(self::Quaternion) option overload;
   return Conjugate(self) . (1/(Norm(self)^2));
  end proc
 ]);

 # QuaternionCommutator computes the commutator of self and x
 export QuaternionCommutator::static := proc(x::Quaternion, y::Quaternion, $)
  return (x . y) - (y . x);
 end proc;

 # display quaternion
 export ModulePrint::static := proc(self::Quaternion, $);
  return cat(self:-real, " + ", self:-i, "i + ", self:-j, "j + ", self:-k, "k"):
 end proc;

 export ModuleApply::static := proc()
  Object(Quaternion, _passed);
 end proc;

 export ModuleCopy::static := proc(new::Quaternion, proto::Quaternion, R::algebraic, imag::algebraic, J::algebraic, K::algebraic, $)
  new:-real := R;
  new:-i := imag;
  new:-j := J;
  new:-k := K;
 end proc;
end module:

q := Quaternion(1, 2, 3, 4):
q1 := Quaternion(2, 3, 4, 5):
q2 := Quaternion(3, 4, 5, 6):
r := 7:

quats := Array([q, q1, q2]):
print("q, q1, q2"):
seq(quats[i], i = 1..3);
print("norms"):
seq(Norm(quats[i]), i = 1..3);
print("negative"):
seq(NegativeQuaternion(quats[i]), i = 1..3);
print("conjugate"):
seq(Conjugate(quats[i]), i = 1..3);
print("addition of real number 7"):
seq(quats[i] + r, i = 1..3);
print("multiplication by real number 7"):
seq(quats[i] . r, i = 1..3);
print("division by real number 7"):
seq(quats[i] / 7, i = 1..3);
print("add quaternions q1 and q2"):
q1 + q2;
print("multiply quaternions q1 and q2");
q1 . q2;
print("multiply quaternions q2 and q1"):
q2 . q1;
print("quaternion commutator of q1 and q2"):
QuaternionCommutator(q1,q2);
print("divide q1 by q2"):
q1 / q2;
Output:

                                 "q, q1, q2"
            1 + 2i + 3j + 4k, 2 + 3i + 4j + 5k, 3 + 4i + 5j + 6k
                                   "norms"
                              1/2     1/2    1/2
                            30   , 3 6   , 86
                                 "negative"
      -1 + -2i + -3j + -4k, -2 + -3i + -4j + -5k, -3 + -4i + -5j + -6k
                                 "conjugate"
        1 + -2i + -3j + -4k, 2 + -3i + -4j + -5k, 3 + -4i + -5j + -6k
                         "addition of real number 7"
            8 + 2i + 3j + 4k, 9 + 3i + 4j + 5k, 10 + 4i + 5j + 6k
                      "multiplication by real number 7"
       7 + 14i + 21j + 28k, 14 + 21i + 28j + 35k, 21 + 28i + 35j + 42k
                         "division by real number 7"
1/7 + 2/7i + 3/7j + 4/7k, 2/7 + 3/7i + 4/7j + 5/7k, 3/7 + 4/7i + 5/7j + 6/7k
                         "add quaternions q1 and q2"
                              5 + 7i + 9j + 11k
                      "multiply quaternions q1 and q2"
                            -56 + 16i + 24j + 26k
                      "multiply quaternions q2 and q1"
                            -56 + 18i + 20j + 28k
                    "quaternion commutator of q1 and q2"
                             0 + -2i + 4j + -2k
                              "divide q1 by q2"
                         34/43 + 1/43i + 0j + 2/43k

Mathematica /Wolfram Language

<<Quaternions`
q=Quaternion[1,2,3,4]
q1=Quaternion[2,3,4,5]
q2=Quaternion[3,4,5,6]
r=7
->Quaternion[1,2,3,4]
->Quaternion[2,3,4,5]
->Quaternion[3,4,5,6]
->7

Abs[q]
->30
-q
->Quaternion[-1,-2,-3,-4]
Conjugate[q]
->Quaternion[1,-2,-3,-4]
r+q
->Quaternion[8,2,3,4]
q+r
->Quaternion[8,2,3,4]
q1+q2
->Quaternion[5,7,9,11]
q*r
->Quaternion[7,14,21,28]
r*q
->Quaternion[7,14,21,28]
q1**q2
->Quaternion[-56,16,24,26]
q2**q1
->Quaternion[-56,18,20,28]

Mercury

A possible implementation of quaternions in Mercury (the simplest representation) would look like this. Note that this is a full module implementation, complete with boilerplate, and that it works by giving an explicit conversion function for floats, converting a float into a quaternion representation of that float. Thus the float value 7.0 gets turned into the quaternion representation q(7.0, 0.0, 0.0, 0.0) through the function call r(7.0).

:- module quaternion.

:- interface.

:- import_module float.

:- type quaternion
    --->    q(  w   :: float,
                i   :: float,
                j   :: float,
                k   :: float    ).

% conversion
:- func r(float) = quaternion is det.

% operations
:- func norm(quaternion) = float is det.
:- func -quaternion = quaternion is det.
:- func conjugate(quaternion) = quaternion is det.
:- func quaternion + quaternion = quaternion is det.
:- func quaternion * quaternion = quaternion is det.

:- implementation.

:- import_module math.

% conversion
r(W) = q(W, 0.0, 0.0, 0.0).

% operations
norm(q(W, I, J, K)) = math.sqrt(W*W + I*I + J*J + K*K).
-q(W, I, J, K) = q(-W, -I, -J, -K).
conjugate(q(W, I, J, K)) = q(W, -I, -J, -K).
q(W0, I0, J0, K0) + q(W1, I1, J1, K1) = q(W0+W1, I0+I1, J0+J1, K0+K1).
q(W0, I0, J0, K0) * q(W1, I1, J1, K1) = q(W0*W1 - I0*I1 - J0*J1 - K0*K1,
                                          W0*I1 + I0*W1 + J0*K1 - K0*J1,
                                          W0*J1 - I0*K1 + J0*W1 + K0*I1,
                                          W0*K1 + I0*J1 - J0*I1 + K0*W1 ).

The following test module puts the module through its paces.

:- module test_quaternion.

:- interface.

:- import_module io.

:- pred main(io::di, io::uo) is det.

:- implementation.

:- import_module quaternion.

:- import_module exception.
:- import_module float.
:- import_module list.
:- import_module string.

:- func to_string(quaternion) = string is det.

main(!IO) :-
    Q  = q(1.0, 2.0, 3.0, 4.0),
    Q1 = q(2.0, 3.0, 4.0, 5.0),
    Q2 = q(3.0, 4.0, 5.0, 6.0),
    R = 7.0,
    QR = r(R),

    io.print("Q = ", !IO), io.print(to_string(Q), !IO), io.nl(!IO),
    io.print("Q1 = ", !IO), io.print(to_string(Q1), !IO), io.nl(!IO),
    io.print("Q2 = ", !IO), io.print(to_string(Q2), !IO), io.nl(!IO),
    io.print("R = ", !IO), io.print(R, !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("1. The norm of a quaternion.\n", !IO),
    io.print("norm(Q) = ", !IO), io.print(norm(Q), !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("2. The negative of a quaternion.\n", !IO),
    io.print("-Q = ", !IO), io.print(to_string(-Q), !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("3. The conjugate of a quaternion.\n", !IO),
    io.print("conjugate(Q) = ", !IO), io.print(to_string(conjugate(Q)), !IO),
        io.nl(!IO),
    io.nl(!IO),

    io.print("4. Addition of a real number and a quaternion.\n", !IO),
    ( Q + QR = QR + Q ->    io.print("Addition is commutative.\n", !IO)
    ;                       io.print("Addition is not commutative.\n", !IO) ),
    io.print("Q + R = ", !IO), io.print(to_string(Q + QR), !IO), io.nl(!IO),
    io.print("R + Q = ", !IO), io.print(to_string(QR + Q), !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("5. Addition of two quaternions.\n", !IO),
    ( Q1 + Q2 = Q2 + Q1 ->  io.print("Addition is commutative.\n", !IO)
    ;                       io.print("Addition is not commutative.\n", !IO) ),
    io.print("Q1 + Q2 = ", !IO), io.print(to_string(Q1 + Q2), !IO), io.nl(!IO),
    io.print("Q2 + Q1 = ", !IO), io.print(to_string(Q2 + Q1), !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("6. Multiplication of a real number and a quaternion.\n", !IO),
    ( Q * QR = QR * Q ->    io.print("Multiplication is commutative.\n", !IO)
    ;                       io.print("Multiplication is not commutative.\n", !IO) ),
    io.print("Q * R = ", !IO), io.print(to_string(Q * QR), !IO), io.nl(!IO),
    io.print("R * Q = ", !IO), io.print(to_string(QR * Q), !IO), io.nl(!IO),
    io.nl(!IO),

    io.print("7. Multiplication of two quaternions.\n", !IO),
    ( Q1 * Q2 = Q2 * Q1 ->  io.print("Multiplication is commutative.\n", !IO)
    ;                       io.print("Multiplication is not commutative.\n", !IO) ),
    io.print("Q1 * Q2 = ", !IO), io.print(to_string(Q1 * Q2), !IO), io.nl(!IO),
    io.print("Q2 * Q1 = ", !IO), io.print(to_string(Q2 * Q1), !IO), io.nl(!IO),
    io.nl(!IO).

to_string(q(I, J, K, W)) = string.format("q(%f, %f, %f, %f)",
                           [f(I), f(J), f(K), f(W)]).
:- end_module test_quaternion.

The output of the above code follows:

% ./test_quaternion 
Q = q(1.000000, 2.000000, 3.000000, 4.000000)
Q1 = q(2.000000, 3.000000, 4.000000, 5.000000)
Q2 = q(3.000000, 4.000000, 5.000000, 6.000000)
R = 7.0

1. The norm of a quaternion.
norm(Q) = 5.477225575051661

2. The negative of a quaternion.
-Q = q(-1.000000, -2.000000, -3.000000, -4.000000)

3. The conjugate of a quaternion.
conjugate(Q) = q(1.000000, -2.000000, -3.000000, -4.000000)

4. Addition of a real number and a quaternion.
Addition is commutative.
Q + R = q(8.000000, 2.000000, 3.000000, 4.000000)
R + Q = q(8.000000, 2.000000, 3.000000, 4.000000)

5. Addition of two quaternions.
Addition is commutative.
Q1 + Q2 = q(5.000000, 7.000000, 9.000000, 11.000000)
Q2 + Q1 = q(5.000000, 7.000000, 9.000000, 11.000000)

6. Multiplication of a real number and a quaternion.
Multiplication is commutative.
Q * R = q(7.000000, 14.000000, 21.000000, 28.000000)
R * Q = q(7.000000, 14.000000, 21.000000, 28.000000)

7. Multiplication of two quaternions.
Multiplication is not commutative.
Q1 * Q2 = q(-56.000000, 16.000000, 24.000000, 26.000000)
Q2 * Q1 = q(-56.000000, 18.000000, 20.000000, 28.000000)

Nim

For simplicity, we have limited the type of quaternion fields to floats (i.e. float64). An implementation could use a generic type in order to allow other field types such as float32.

import math, tables

type Quaternion* = object
  a, b, c, d: float

func initQuaternion*(a, b, c, d = 0.0): Quaternion =
  Quaternion(a: a, b: b, c: c, d: d)

func `-`*(q: Quaternion): Quaternion =
  initQuaternion(-q.a, -q.b, -q.c, -q.d)

func `+`*(q: Quaternion; r: float): Quaternion =
  initQuaternion(q.a + r, q.b, q.c, q.d)

func `+`*(r: float; q: Quaternion): Quaternion =
  initQuaternion(q.a + r, q.b, q.c, q.d)

func `+`*(q1, q2: Quaternion): Quaternion =
  initQuaternion(q1.a + q2.a, q1.b + q2.b, q1.c + q2.c, q1.d + q2.d)

func `*`*(q: Quaternion; r: float): Quaternion =
  initQuaternion(q.a * r, q.b * r, q.c * r, q.d * r)

func `*`*(r: float; q: Quaternion): Quaternion =
  initQuaternion(q.a * r, q.b * r, q.c * r, q.d * r)

func `*`*(q1, q2: Quaternion): Quaternion =
  initQuaternion(q1.a * q2.a - q1.b * q2.b - q1.c * q2.c - q1.d * q2.d,
                 q1.a * q2.b + q1.b * q2.a + q1.c * q2.d - q1.d * q2.c,
                 q1.a * q2.c - q1.b * q2.d + q1.c * q2.a + q1.d * q2.b,
                 q1.a * q2.d + q1.b * q2.c - q1.c * q2.b + q1.d * q2.a)

func conjugate*(q: Quaternion): Quaternion =
  initQuaternion(q.a, -q.b, -q.c, -q.d)

func norm*(q: Quaternion): float =
  sqrt(q.a * q.a + q.b * q.b + q.c * q.c + q.d * q.d)

func `==`*(q: Quaternion; r: float): bool =
  if q.b != 0 or q.c != 0 or q.d != 0: false
  else: q.a == r

func `$`(q: Quaternion): string =
  ## Return the representation of a quaternion.
  const Letter = {"a": "", "b": "i", "c": "j", "d": "k"}.toTable
  if q == 0: return "0"
  for name, value in q.fieldPairs:
    if value != 0:
      var val = value
      if result.len != 0:
        result.add if value >= 0: '+' else: '-'
        val = abs(val)
      result.add $val & Letter[name]


when isMainModule:
  let
    q = initQuaternion(1, 2, 3, 4)
    q1 = initQuaternion(2, 3, 4, 5)
    q2 = initQuaternion(3, 4, 5, 6)
    r = 7.0

  echo "∥q∥ = ", norm(q)
  echo "-q = ", -q
  echo "q* = ", conjugate(q)
  echo "q + r = ", q + r
  echo "r + q = ", r + q
  echo "q1 + q2 = ", q1 + q2
  echo "qr = ", q * r
  echo "rq = ", r * q
  echo "q1 * q2 = ", q1 * q2
  echo "q2 * q1 = ", q2 * q1
Output:
∥q∥ = 5.477225575051661
-q = -1.0-2.0i-3.0j-4.0k
q* = 1.0-2.0i-3.0j-4.0k
q + r = 8.0+2.0i+3.0j+4.0k
r + q = 8.0+2.0i+3.0j+4.0k
q1 + q2 = 5.0+7.0i+9.0j+11.0k
qr = 7.0+14.0i+21.0j+28.0k
rq = 7.0+14.0i+21.0j+28.0k
q1 * q2 = -56.0+16.0i+24.0j+26.0k
q2 * q1 = -56.0+18.0i+20.0j+28.0k

As can be seen, q1 * q2 != q2 * q1.

OCaml

This implementation was build strictly to the specs without looking (too much) at other implementations. The implementation as a record type with only floats is said (on the ocaml mailing list) to be especially efficient. Put this into a file quaternion.ml:

type quaternion = {a: float; b: float; c: float; d: float}

let norm q = sqrt (q.a**2.0 +.
                   q.b**2.0 +.
                   q.c**2.0 +.
                   q.d**2.0 )

let floatneg r = ~-. r  (* readability *)

let negative q =
  {a = floatneg q.a; 
   b = floatneg q.b;
   c = floatneg q.c;
   d = floatneg q.d }

let conjugate q =
  {a = q.a; 
   b = floatneg q.b;
   c = floatneg q.c;
   d = floatneg q.d }

let addrq r q = {q with a = q.a +. r}

let addq q1 q2 =
  {a = q1.a +. q2.a;
   b = q1.b +. q2.b;
   c = q1.c +. q2.c;
   d = q1.d +. q2.d  }

let multrq r q =
  {a = q.a *. r;
   b = q.b *. r;
   c = q.c *. r;
   d = q.d *. r  }
   
let multq q1 q2 =  
        {a = q1.a*.q2.a -. q1.b*.q2.b -. q1.c*.q2.c -. q1.d*.q2.d;
         b = q1.a*.q2.b +. q1.b*.q2.a +. q1.c*.q2.d -. q1.d*.q2.c;
         c = q1.a*.q2.c -. q1.b*.q2.d +. q1.c*.q2.a +. q1.d*.q2.b;
         d = q1.a*.q2.d +. q1.b*.q2.c -. q1.c*.q2.b +. q1.d*.q2.a  }
   
let qmake a b c d = {a;b;c;d} (* readability omitting a= b=... *)

let qstring q =
  Printf.sprintf "(%g, %g, %g, %g)" q.a q.b q.c q.d ;;

(* test data *)  
let q  = qmake 1.0  2.0  3.0  4.0
let q1 = qmake 2.0  3.0  4.0  5.0
let q2 = qmake 3.0  4.0  5.0  6.0
let r  = 7.0

let () = (* written strictly to spec *)
  let pf = Printf.printf in
  pf "starting with data q=%s, q1=%s,  q2=%s, r=%g\n" (qstring q) (qstring q1) (qstring q2) r;
  pf "1. norm of      q     = %g \n" (norm q) ;
  pf "2. negative of  q     = %s \n" (qstring (negative q));
  pf "3. conjugate of q     = %s \n" (qstring (conjugate q));
  pf "4. adding r to q      = %s \n" (qstring (addrq r q));
  pf "5. adding q1 and q2   = %s \n" (qstring (addq q1 q2));
  pf "6. multiply r and q   = %s \n" (qstring (multrq r q));
  pf "7. multiply q1 and q2 = %s \n" (qstring (multq q1 q2));
  pf "8. instead q2 * q1    = %s \n" (qstring (multq q2 q1));
  pf "\n";

using this file on the command line will produce:

$ ocaml quaternion.ml 
starting with data q=(1, 2, 3, 4), q1=(2, 3, 4, 5),  q2=(3, 4, 5, 6), r=7
1. norm of      q     = 5.47723 
2. negative of  q     = (-1, -2, -3, -4) 
3. conjugate of q     = (1, -2, -3, -4) 
4. adding r to q      = (8, 2, 3, 4) 
5. adding q1 and q2   = (5, 7, 9, 11) 
6. multiply r and q   = (7, 14, 21, 28) 
7. multiply q1 and q2 = (-56, 16, 24, 26) 
8. instead q2 * q1    = (-56, 18, 20, 28) 

For completeness, and since data types are of utmost importance in OCaml, here the types produced by pasting the code into the toplevel (ocaml is the toplevel):

type quaternion = { a : float; b : float; c : float; d : float; }
val norm : quaternion -> float = <fun>
val floatneg : float -> float = <fun>
val negative : quaternion -> quaternion = <fun>
val conjugate : quaternion -> quaternion = <fun>
val addrq : float -> quaternion -> quaternion = <fun>
val addq : quaternion -> quaternion -> quaternion = <fun>
val multrq : float -> quaternion -> quaternion = <fun>
val multq : quaternion -> quaternion -> quaternion = <fun>
val qmake : float -> float -> float -> float -> quaternion = <fun>
val qstring : quaternion -> string = <fun>

Octave

There is an add-on package (toolbox) to Octave available from http://octave.sourceforge.net/quaternion/

Such a package can be install with the command:

pkg install -forge quaternion

Here is a sample interactive session solving the task:

> q = quaternion (1, 2, 3, 4)
q = 1 + 2i + 3j + 4k
>  q1 = quaternion (2, 3, 4, 5)
q1 = 2 + 3i + 4j + 5k
> q2 = quaternion (3, 4, 5, 6)
q2 = 3 + 4i + 5j + 6k
> r = 7
r =  7
> norm(q)
ans =  5.4772
> -q
ans = -1 - 2i - 3j - 4k
> conj(q)
ans = 1 - 2i - 3j - 4k
> q + r
ans = 8 + 2i + 3j + 4k
> q1 + q2
ans = 5 + 7i + 9j + 11k
> q * r
ans = 7 + 14i + 21j + 28k
> q1 * q2
ans = -56 + 16i + 24j + 26k
> q1 == q2
ans = 0

Oforth

Setting a priority (here 160) to Quaternion class and defining #asQuaternion, integers and floats can be fully mixed with quaternions. neg is defined as "0 self -" into Number class, so no need to define it (if #- is defined).

160 Number Class newPriority: Quaternion(a, b, c, d)

Quaternion method: _a  @a ;
Quaternion method: _b  @b ;
Quaternion method: _c  @c ;
Quaternion method: _d  @d ;

Quaternion method: initialize  := d := c := b := a ;
Quaternion method: <<  '(' <<c @a << ',' <<c @b << ',' <<c @c << ',' <<c @d << ')' <<c ;

Integer method: asQuaternion  self 0 0 0 Quaternion new ;
Float   method: asQuaternion  self 0 0 0 Quaternion new ;

Quaternion method: ==(q)  q _a @a == q _b @b == and q _c @c == and q _d @d == and ;
Quaternion method: norm   @a sq @b sq + @c sq + @d sq + sqrt ;
Quaternion method: conj   @a  @b neg  @c neg  @d neg Quaternion new ;
Quaternion method: +(q)   Quaternion new(q _a @a +, q _b @b +, q _c @c +, q _d @d +) ;
Quaternion method: -(q)   Quaternion new(q _a @a -, q _b @b -, q _c @c -, q _d @d -) ;

Quaternion method: *(q) 
   Quaternion new(q _a @a * q _b @b * - q _c @c * - q _d @d * -,
                  q _a @b * q _b @a * + q _c @d * + q _d @c * -,
                  q _a @c * q _b @d * - q _c @a * + q _d @b * +,
                  q _a @d * q _b @c * + q _c @b * - q _d @a * + ) ;

Usage :

: test
| q q1 q2 r |

   Quaternion new(1, 2, 3, 4) ->q
   Quaternion new(2, 3, 4, 5) ->q1
   Quaternion new(3, 4, 5, 6) ->q2
   7.0 -> r

   System.Out "q       = " << q << cr
   System.Out "q1      = " << q1 << cr
   System.Out "q2      = " << q2 << cr

   System.Out "norm q  = " << q norm << cr 
   System.Out "neg q   = " << q neg << cr 
   System.Out "conj q  = " << q conj << cr
   System.Out "q +r    = " << q r + << cr
   System.Out "q1 + q2 = " << q1 q2 + << cr
   System.Out "q * r   = " << q r * << cr
   System.Out "q1 * q2 = " << q1 q2 * << cr
   q1 q2 * q2 q1 * == ifFalse: [ "q1q2 and q2q1 are different quaternions" println ] ;
Output:
q       = (1,2,3,4)
q1      = (2,3,4,5)
q2      = (3,4,5,6)
norm q  = 5.47722557505166
neg q   = (-1,-2,-3,-4)
conj q  = (1,-2,-3,-4)
q +r    = (8,2,3,4)
q1 + q2 = (5,7,9,11)
q * r   = (7,14,21,28)
q1 * q2 = (-56,16,24,26)
q1q2 and q2q1 are different quaternions

Ol

See also the entry for Scheme.

;;
;; This program is written to run without modification both in Otus
;; Lisp and in any of many Scheme dialects. I assume the presence of
;; "case-lambda", but not of "let-values". The program has worked
;; (without modification) in Otus Lisp 2.4, Guile >= 2.0 (but not in
;; Guile version 1.8), CHICKEN Scheme 5.3.0, Chez Scheme 9.5.8, Gauche
;; Scheme 0.9.12, Ypsilon 0.9.6-update3.
;;
;; Here a quaternion is represented as a linked list of four real
;; numbers. Such a representation probably has the greatest
;; portability between Scheme dialects. However, this representation
;; can be replaced, simply by redefining the procedures "quaternion?",
;; "quaternion-components", "quaternion->list", and "quaternion".
;;

(define (quaternion? q)               ; Can q be used as a quaternion?
  (and (pair? q)
       (let ((a (car q))
             (q (cdr q)))
         (and (real? a) (pair? q)
              (let ((b (car q))
                    (q (cdr q)))
                (and (real? b) (pair? q)
                     (let ((c (car q))
                           (q (cdr q)))
                       (and (real? c) (pair? q)
                            (let ((d (car q))
                                  (q (cdr q)))
                              (and (real? d) (null? q)))))))))))

(define (quaternion-components q)       ; Extract the basis components.
  (let ((a (car q))
        (q (cdr q)))
    (let ((b (car q))
          (q (cdr q)))
      (let ((c (car q))
            (q (cdr q)))
        (let ((d (car q)))
          (values a b c d))))))

(define (quaternion->list q)     ; Get a list of the basis components.
  q)

(define quaternion                      ; Make a quaternion.
  (case-lambda
    ((a b c d)
     ;; Make the quaternion from basis components.
     (list a b c d))
    ((q)
     ;; Make the quaternion from a scalar or from another quaternion.
     ;; WARNING: in the latter case, the quaternion is NOT
     ;; copied. This is not a problem, if you avoid things like
     ;; "set-car!" and "set-cdr!".
     (if (real? q)
         (list q 0 0 0)
         q))))

(define (quaternion-norm q)      ; The euclidean norm of a quaternion.
  (let ((q (quaternion q)))
    (call-with-values (lambda () (quaternion-components q))
      (lambda (a b c d)
        (sqrt (+ (* a a) (* b b) (* c c) (* d d)))))))

(define (quaternion-conjugate q)        ; Conjugate a quaternion.
  (let ((q (quaternion q)))
    (call-with-values (lambda () (quaternion-components q))
      (lambda (a b c d)
        (quaternion a (- b) (- c) (- d))))))

(define quaternion+                     ; Add quaternions.
  (let ((quaternion-add
         (lambda (q1 q2)
           (let ((q1 (quaternion q1))
                 (q2 (quaternion q2)))
             (call-with-values
                 (lambda () (quaternion-components q1))
               (lambda (a1 b1 c1 d1)
                 (call-with-values
                     (lambda () (quaternion-components q2))
                   (lambda (a2 b2 c2 d2)
                     (quaternion (+ a1 a2) (+ b1 b2)
                                 (+ c1 c2) (+ d1 d2))))))))))
    (case-lambda
      (() (quaternion 0))