Jump to content

Quaternion type

From Rosetta Code
Revision as of 20:53, 27 November 2010 by rosettacode>Mr2001 (added C# solution)
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 written sometimes 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 this numbering system, 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 Description
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.

Your task is to 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 your language has built-in support for quaternions then use it.

Ada

The package specification (works with any floating-point type): <lang Ada>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;</lang> The package implementation: <lang Ada>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;</lang> Test program: <lang Ada>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;</lang> Sample 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 - no extensions to language used
Works with: ALGOL 68G version Any - tested with release 1.18.0-9h.tiny

<lang algol68>MODE QUAT = STRUCT(REAL re, i, j, k); MODE QUATERNION = QUAT; MODE SUBQUAT = UNION(QUAT, #COMPL, # REAL#, INT, [4]REAL, [4]INT # );

MODE CLASSQUAT = STRUCT( 

   PROC (REF QUAT #new#, REAL #re#, REAL #i#, REAL #j#, REAL #k#)REF QUAT new,
   PROC (REF QUAT #self#)QUAT conjugate,
   PROC (REF QUAT #self#)REAL norm sq,
   PROC (REF QUAT #self#)REAL 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, REAL re, i, j, k)REF QUAT: (
       # 'Defaults all parts of quaternion to zero' #
       IF new ISNT REF QUAT(NIL) THEN new ELSE HEAP QUAT FI := (re, i, j, k)
   ),

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

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

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

 # PROC reciprocal =#(REF QUAT self)QUAT:(
       REAL n2 = (norm sq OF class quat)(self);
       QUAT conj = (conjugate OF class quat)(self);
       (re 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, re OF self>=0, re 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:
       (-re OF self, -i OF self, -j OF self, -k OF self),

 # PROC add =#(REF QUAT self, SUBQUAT other)QUAT:
       CASE other IN
           (QUAT other): (re OF self + re OF other, i OF self + i OF other, j OF self + j OF other, k OF self + k OF other),
           (REAL other): (re OF self + other, i OF self, j OF self, k OF self)
       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): (re OF self - re OF other, i OF self - i OF other, j OF self - j OF other, k OF self - k OF other),
           (REAL other): (re OF self - other, i OF self, j OF self, k OF self)
       ESAC,

 # PROC mul =#(REF QUAT self, SUBQUAT other)QUAT:
       CASE other IN
           (QUAT other):(
                re OF self*re OF other - i OF self*i  OF other - j OF self*j  OF other - k OF self*k  OF other,
                re OF self*i  OF other + i OF self*re OF other + j OF self*k  OF other - k OF self*j  OF other,
                re OF self*j  OF other - i OF self*k  OF other + j OF self*re OF other + k OF self*i  OF other,
                re OF self*k  OF other + i OF self*j  OF other - j OF self*i  OF other + k OF self*re OF other
           ),
           (REAL other): ( re OF self * other, i OF self * other, j OF self * other, k OF self * other)
       ESAC,

 # PROC rmul =#(REF QUAT self, SUBQUAT other)QUAT:
       CASE other IN
         (QUAT other): (mul OF class quat)(LOC QUAT := other, self),
         (REAL other): (mul OF class quat)(self, other)
       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)),
           (REAL other): (mul OF class quat)(self, 1/other)
       ESAC,

 # PROC rdiv =#(REF QUAT self, SUBQUAT other)QUAT:
       CASE other IN
         (QUAT other): (div OF class quat)(LOC QUAT := other, self),
         (REAL other): (div OF class quat)(LOC QUAT := (other, 0, 0, 0), self)
       ESAC,

 # PROC exp =#(REF QUAT self)QUAT: (
   QUAT fac := self;
   QUAT sum := 1.0 + fac;
   FOR i FROM 2 WHILE ABS(fac + small real) /= small real DO
     VOID(sum +:= (fac *:= self / REAL(i)))
   OD;
   sum
 )

);

FORMAT real fmt = $g(-0, 4)$; FORMAT signed fmt = $b("+", "")f(real fmt)$;

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

PRIO INIT = 1; OP INIT = (REF QUAT new)REF QUAT: new := (0, 0, 0, 0); OP INIT = (REF QUAT new, []REAL rijk)REF QUAT: 

   (new OF class quat)(LOC QUAT := new, rijk[1], rijk[2], rijk[3], rijk[4]);

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)REAL:   (norm OF class quat)(LOC QUAT := q),
  REPR = (QUAT q)STRING: (repr OF class quat)(LOC QUAT := q);
  1. missing: Diadic: I, J, K END #

OP +:= = (REF QUAT a, QUAT b)QUAT: a:=( add OF class quat)(a, b), 

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

OP -:= = (REF QUAT a, QUAT b)QUAT: a:=( sub OF class quat)(a, b), 

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

PRIO *=: = 1, /=: = 1; OP *:= = (REF QUAT a, QUAT b)QUAT: a:=( mul OF class quat)(a, b), 

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

OP /:= = (REF QUAT a, QUAT b)QUAT: a:=( div OF class quat)(a, b), 

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

OP + = (QUAT a, b)QUAT: ( add OF class quat)(LOC QUAT := a, b), 

  + = (QUAT a, REAL b)QUAT: ( add OF class quat)(LOC QUAT := a, b),
  + = (REAL 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, REAL b)QUAT: ( sub OF class quat)(LOC QUAT := a, b),
  - = (REAL 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, REAL b)QUAT: ( mul OF class quat)(LOC QUAT := a, b),
  * = (REAL 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, REAL b)QUAT: ( div OF class quat)(LOC QUAT := a, b),
  / = (REAL a, QUAT b)QUAT: ( div OF class quat)(LOC QUAT := b, 1/a);

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

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

test:( 

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

   printf((
       $"r = "      f(real fmt)l$, r,
       $"q = "      f(quat fmt)l$, q,
       $"q1 = "     f(quat fmt)l$, q1,
       $"q2 = "     f(quat fmt)l$, q2,
       $"ABS q = "  f(real fmt)", "$, ABS q,
       $"ABS q1 = " f(real fmt)", "$, ABS q1,
       $"ABS q2 = " f(real 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 

   QUAT i=(0, 1, 0, 0),
        j=(0, 0, 1, 0),
        k=(0, 0, 0, 1);

   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 * i) = " f(quat fmt)l$, quat exp(pi * k)
   ));
   print((REPR(-q1*q2), ", ", REPR(-q2*q1), new line))

)</lang> 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 = -46.0000+12.0000i+16.0000j+20.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 * i) = -1.0000+.0000i+.0000j+.0000k
+56.0000-16.0000i-24.0000j-26.0000k, +56.0000-18.0000i-20.0000j-28.0000k

C

<lang c>#include <stdio.h>

  1. include <stdlib.h>
  2. include <stdbool.h>
  3. 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;

}


  1. define A(N) (a->q[(N)])
  2. define B(N) (b->q[(N)])
  3. 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);

}

  1. undef A
  2. undef B
  3. 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]); }</lang>

<lang c>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;

}</lang>

C#

<lang csharp>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

}</lang>

Demonstration: <lang csharp>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) ? "==" : "!=");
   }

}</lang>

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

D

Works with: D version 2.049

<lang d>module quaternion ; import std.complex ; import std.math ; import std.numeric; import std.traits; import std.string ;

struct H(T) if (isFloatingPoint!T) { 

   alias Complex!(T) CT ;
   union {
       struct { T re, i, j, k ;} // default init to NaN
       struct { CT x, y ; }
       struct { T[4] vector ; }
   }

   string toString() { return format("%s", vector) ; }

   @property T norm2() const {    // norm squared
       return re*re + i*i + j*j + k*k ;
   }
   @property T abs() const {      // norm
       return sqrt(norm2) ;
   }
   @property T arg() const {      // theta -- this may be incorrect...
       return acos(re / abs) ;
   }
   @property H!T conj() const {   // conjugate
       return H!T(re, -i, -j, -k);
   }
   @property H!T recip() const {  // reciprocal
       return H!T(re/norm2, -i/norm2, -j/norm2, -k/norm2);
   }
   @property H!T pureim() const { // pure imagery
       return H!T(0, i, j, k);
   }
   @property H!T versor() const { // unit versor
       return this/abs;
   }
   @property H!T iversor() const {// unit versor of imagery part
       return pureim/pureim.abs;
   }
   // assignment
   H!T opAssign(U : T)(H!U z) {
       x = z.x ;  y = z.y ;
       return this;
   }
   H!T opAssign(U : T)(Complex!U c) {
       x = c ;  y = 0 ;
       return this;
   }
   H!T opAssign(U : T)(U r) if (isNumeric!U) {
       re = r ; i = 0 ; y = 0 ;
       return this;
   }
   // test for equal, not ordered so no opCmp
   bool opEquals(U : T)(H!U z) const {
       return re == z.re && i == z.i && j == z.j && k == z.k ;
   }
   bool opEquals(U : T)(Complex!U c) const {
       return re == c.re && i == c.im && j == 0 && k == 0 ;
   }
   bool opEquals(U : T)(U r) const if (isNumeric!U) {
       return re == r && i == 0 && j == 0 && k == 0 ;
   }
   // unary op
   H!T opUnary(string op)() const if (op == "+") {
       return this;
   }
   H!T opUnary(string op)() const if (op == "-") {
       return H!T(-re, -i, -j, -k);
   }
   // binary op, Quaternion on left of op
   H!(CommonType!(T,U)) opBinary(string op, U)(H!U z) const  {
       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
   H!(CommonType!(T,U)) opBinary(string op, U)(Complex!U c) const  {
       alias typeof(return) C;
       return opBinary!(op)(C(c.re, c.im, 0, 0));
   }
   // for scalar
   H!(CommonType!(T,U)) opBinary(string op, U)(U r) const
       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) ;
       }
       static if(op == "^^") {
           return pow(r) ;
       }
   }
   // power function
   H!(CommonType!(T,U)) pow(U)(U r) const if (isNumeric!U) {
       return (abs^^r) * exp(r * iversor * arg) ;
   }
   // handle binary op if Quaternion on right of op and left is not quaternion
   H!(CommonType!(T,U)) opBinaryRight(string op, U)(Complex!U c)
       const {
       alias typeof(return) C;
       auto w = C(c.re, c.im, 0, 0) ;
       return w.opBinary!(op)(this);
   }
   H!(CommonType!(T,U)) opBinaryRight(string op, U)(U r) const
       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 ;
       }
   }

} // exp / log functions HT exp(HT)(HT z) if (is(HT T == H!T)) { 

   auto inorm = z.pureim.abs ;
   return std.math.exp(z.re)*(cos(inorm) + z.iversor * sin(inorm)) ;

} HT log(HT)(HT z) if (is(HT T == H!T)) { 

   return std.math.log(z.abs) + z.iversor*acos(z.re/z.abs) ;

}

// test program

import std.stdio ; void main() { 

   H!real q, q1, q2 ;
   real r = 7 ;
   q.vector = [1,2,3,4] ;
   q1 = H!real(2,3,4,5) ;
   q2 = H!real(3,4,5,6) ;

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

   H!real i,j,k ;
   i.vector = [0,1,0,0] ; j.vector = [0,0,1,0] ; k.vector = [0,0,0,1] ;
   writefln("9.1 i*i, J*j, k*k, i*j*k : %s, %s, %s, %s",
       i*i, j*j, k*k, i*j*k) ;
   writefln("9.2              q1 / q2 : %s", q1 / q2) ;
   writefln("9.3         q1 / q2 * q2 : %s", q1 / q2 * q2) ;
   writefln("            q2 * q1 / q2 : %s", q2 * q1 / q2) ;
   writefln("9.4          exp(pi * i) : %s", exp(PI*i)) ;
   writefln("             exp(pi * j) : %s", exp(PI*j)) ;
   writefln("             exp(pi * k) : %s", exp(PI*k)) ;
   writefln("                  exp(q) : %s", exp(q)) ;
   writefln("                  log(q) : %s", log(q)) ;
   writefln("             exp(log(q)) : %s", exp(log(q))) ;
   writefln("             log(exp(q)) : %s", log(exp(q))) ;
   auto s = log(exp(q)) ;
   writefln("9.5  let s = log(exp(q)) : %s", s) ;
   writefln("                  exp(s) : %s", exp(s)) ;
   writefln("                  log(s) : %s", log(s)) ;
   writefln("             exp(log(s)) : %s", exp(log(s))) ;
   writefln("             log(exp(s)) : %s", log(exp(s))) ;

} /+ outputs 1. q - norm  : 5.47723 2. q - negative  : [-1,-2,-3,-4] 3. q - conjugate : [1,-2,-3,-4] 4. r + q  : [8,2,3,4] 

  q + r         : [8,2,3,4]

5. q1 + q2  : [5,7,9,11] 6. r * q  : [7,14,21,28] 

  q * r         : [7,14,21,28]

7. q1 * q2  : [-56,16,24,26] 

  q2 * q1       : [-56,18,20,28]

8. q1q2 != q2Q1  ? true 9.1 i*i, J*j, k*k, i*j*k : [-1,0,0,0], [-1,0,0,0], [-1,0,0,0], [-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) : [1.69392,-0.78956,-1.18434,-1.57912]
                 log(q) : [1.7006,0.51519,0.772785,1.03038]
            exp(log(q)) : [1,2,3,4]
            log(exp(q)) : [1,-0.333516,-0.500275,-0.667033]

9.5 let s = log(exp(q)) : [1,-0.333516,-0.500275,-0.667033] 

                 exp(s) : [1.69392,-0.78956,-1.18434,-1.57912]
                 log(s) : [0.295679,-0.271754,-0.407631,-0.543508]
            exp(log(s)) : [1,-0.333516,-0.500275,-0.667033]
            log(exp(s)) : [1,-0.333516,-0.500275,-0.667033]

+/</lang>

Forth

<lang 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)</lang>

Fortran

Works with: Fortran version 90 and later

<lang fortran>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</lang> 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

Haskell

<lang haskell>import Control.Monad import Control.Arrow import Data.List

data Quaternion = Q Double Double Double Double 

 deriving (Show, Ord, Eq)

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

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

quaternionFromScalar s = Q s 0 0 0

listFromQ (Q a b c d) = [a,b,c,d] quaternionFromList [a, b, c, d] = Q a b c d

addQ, subQ, mulQ :: Quaternion -> Quaternion -> Quaternion addQ (Q a b c d) (Q p q r s) = Q (a+p) (b+q) (c+r) (d+s)

subQ (Q a b c d) (Q p q r s) = Q (a-p) (b-q) (c-r) (d-s)

mulQ (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)

normQ = sqrt. sum. join (zipWith (*)). listFromQ

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

negQ (Q a b c d) = Q (-a) (-b) (-c) (-d)</lang> To use with the Examples: <lang haskell>[q,q1,q2] = map quaternionFromList [[1..4],[2..5],[3..6]] -- a*b == b*a test :: Quaternion -> Quaternion -> Bool test a b = a `mulQ` b == b `mulQ` a</lang> Examples: 

*Main> mulQ (Q 0 1 0 0) $ mulQ (Q 0 0 1 0) (Q 0 0 0 1) -- i*j*k
Q (-1.0) 0.0 0.0 0.0

*Main> test q1 q2
False

*Main> mulQ q1 q2
Q (-56.0) 16.0 24.0 26.0

*Main> flip mulQ q1 q2
Q (-56.0) 18.0 20.0 28.0

*Main> imagQ q
[2.0,3.0,4.0]

J

Derived from the j wiki:

<lang j> NB. utilities 

  ip=:   +/ .*             NB. inner product
  t4=. (_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 t4 ip ])&toQ  NB. x * y</lang>

Example use:

<lang> 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</lang>

Lua

<lang 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</lang> Examples: <lang lua>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</lang>

OCaml

The implementation as a file q.ml: <lang ocaml>type quaternion = float * float * float * float

let q a b c d = (a, b, c, d)

let to_real (r, _, _, _) = r let imag (_, i, j, k) = (i, j, k)

let quaternion_of_scalar s = (s, 0.0, 0.0, 0.0)

let to_list (a, b, c, d) = [a; b; c; d] let of_list = function [a; b; c; d] -> (a, b, c, d) 

 | _ -> invalid_arg "of_list"

let ( + ) = ( +. ) let ( - ) = ( -. ) let ( * ) = ( *. ) let ( / ) = ( /. )

let addr (a, b, c, d) r = (a+r, b, c, d) let mulr (a, b, c, d) r = (a*r, b*r, c*r, d*r)

let add (a, b, c, d) (p, q, r, s) = (a+p, b+q, c+r, d+s)

let sub (a, b, c, d) (p, q, r, s) = (a-p, b-q, c-r, d-s)

let mul (a, b, c, d) (p, q, r, s) = 

 ( 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 )

let norm2 (a, b, c, d) = 

 ( a * a +
   b * b +
   c * c +
   d * d )

let norm q = sqrt(norm2 q)

let conj (a, b, c, d) = (a, -. b, -. c, -. d) let neg (a, b, c, d) = (-. a, -. b, -. c, -. d)

let unit ((a, b, c, d) as q) = 

 let n = norm q in
 (a/n, b/n, c/n, d/n)

let reciprocal ((a, b, c, d) as q) = 

 let n2 = norm2 q in
 (a/n2, b/n2, c/n2, d/n2)</lang>

and the interface as a file q.mli: <lang ocaml>type quaternion = float * float * float * float

val q : float -> float -> float -> float -> quaternion val to_real : quaternion -> float val imag : quaternion -> float * float * float val quaternion_of_scalar : float -> quaternion val to_list : quaternion -> float list val of_list : float list -> quaternion val addr : quaternion -> float -> quaternion val mulr : quaternion -> float -> quaternion val add : quaternion -> quaternion -> quaternion val sub : quaternion -> quaternion -> quaternion val mul : quaternion -> quaternion -> quaternion val norm : quaternion -> float val conj : quaternion -> quaternion val neg : quaternion -> quaternion val unit : quaternion -> quaternion val reciprocal : quaternion -> quaternion</lang>

using this module in the interactive interpreter: 

$ ocamlc -c q.mli 
$ ocamlc -c q.ml 
$ ocaml q.cmo
        Objective Caml version 3.11.2

# open Q ;;
# let q1 = q 2.0 3.0 4.0 5.0
  and q2 = q 3.0 4.0 5.0 6.0 ;;
val q1 : Q.quaternion = (2., 3., 4., 5.)
val q2 : Q.quaternion = (3., 4., 5., 6.)
# (mul q1 q2) <> (mul q2 q1) ;;
- : bool = true

Perl 6

<lang perl6>class Quaternion { 

   has Real ( $.r, $.i, $.j, $.k );

   multi method new ( Real $r, Real $i, Real $j, Real $k ) {
       self.bless: *, :$r, :$i, :$j, :$k;
   }
   method Str   (  ) { "$.r + {$.i}i + {$.j}j + {$.k}k" }
   method reals (  ) { $.r, $.i, $.j, $.k }
   method conj  (  ) { self.new: $.r, -$.i, -$.j, -$.k }
   method norm  (  ) { sqrt [+] self.reals X** 2 }

}

subset Qu of Quaternion; # Makes a short alias

multi sub infix:<eqv> ( Qu $a, Qu $b ) { [and] $a.reals Z== $b.reals } multi sub infix:<+> ( Qu $a, Real $b ) { $a.new: $b+$a.r, $a.i, $a.j, $a.k } multi sub infix:<+> ( Real $a, Qu $b ) { $b.new: $a+$b.r, $b.i, $b.j, $b.k } multi sub infix:<+> ( Qu $a, Qu $b ) { $a.new: |( $a.reals Z+ $b.reals ) } multi sub prefix:<-> ( Qu $a ) { $a.new: |( $a.reals X* -1 ) } multi sub infix:<*> ( Qu $a, Real $b ) { $a.new: |( $a.reals X* $b ) } multi sub infix:<*> ( Real $a, Qu $b ) { $b.new: |( $b.reals X* $a ) } multi sub infix:<*> ( Qu $a, Qu $b ) { 

   my @a_rijk            = $a.reals;
   my ( $r, $i, $j, $k ) = $b.reals;
   return $a.new: ( [+] @a_rijk Z* $r, -$i, -$j, -$k ), # real
                  ( [+] @a_rijk Z* $i,  $r,  $k, -$j ), # i
                  ( [+] @a_rijk Z* $j, -$k,  $r,  $i ), # j
                  ( [+] @a_rijk Z* $k,  $j, -$i,  $r ); # k

}

my Quaternion $q .= new: 1, 2, 3, 4; my Quaternion $q1 .= new: 2, 3, 4, 5; my Quaternion $q2 .= new: 3, 4, 5, 6; my $r = 7;

say "1) q norm = {$q.norm}"; say "2) -q = {-$q}"; say "3) q conj = {$q.conj}"; say "4) q + r = {$q + $r}"; say "5) q1 + q2 = {$q1 + $q2}"; say "6) q * r = {$q * $r}"; say "7) q1 * q2 = {$q1 * $q2}"; say "8) q1q2 { $q1 * $q2 eqv $q2 * $q1 ?? '==' !! '!=' } q2q1"; </lang>

Output:

1) q norm  = 5.47722557505166
2) -q      = -1 + -2i + -3j + -4k
3) q conj  = 1 + -2i + -3j + -4k
4) q  + r  = 8 + 2i + 3j + 4k
5) q1 + q2 = 5 + 7i + 9j + 11k
6) q  * r  = 7 + 14i + 21j + 28k
7) q1 * q2 = -56 + 16i + 24j + 26k
8) q1q2 != q2q1

PureBasic

<lang PureBasic>Structure Quaternion 

 a.f: b.f
 c.f: d.f

EndStructure

Structure Quaternion2 

 Qa.Quaternion: Qb.Quaternion
 Qc.Quaternion: Qd.Quaternion

EndStructure

Procedure.f QNorm(*Q.Quaternion) 

 Protected Result.f
 Result=Sqr(Pow(*Q\a,2)+ Pow(*Q\b,2)+ Pow(*Q\c,2)+ Pow(*Q\d,2))
 ProcedureReturn Result

EndProcedure

Procedure.i QNeg(*Q.Quaternion) 

 Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion))
 If *B
   With *Q
     *B\a=-\a: *B\b=-\b
     *B\c=-\c: *B\d=-\d
   EndWith
 EndIf
 ProcedureReturn *B

EndProcedure

Procedure.i QConj(*Q.Quaternion) 

 Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion))
 If *B
   With *Q
     *B\a= \a: *B\b=-\b
     *B\c=-\c: *B\d=-\d
   EndWith
 EndIf
 ProcedureReturn *B

EndProcedure

Procedure.i QAddReal(R.f, *Q.Quaternion) 

 Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion))
 If *B
   With *Q
     *B\a=R+\a: *B\b= -\b
     *B\c= -\c: *B\d= -\d
   EndWith
 EndIf
 ProcedureReturn *B

EndProcedure

Procedure.i QAddQuaternion(*Q1.Quaternion, *Q2.Quaternion) 

 Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion))
 If *B
   *B\a=*Q1\a + *Q2\a: *B\b=*Q1\b + *Q2\b
   *B\c=*Q1\c + *Q2\c: *B\d=*Q1\d + *Q2\d
 EndIf
 ProcedureReturn *B

EndProcedure

Procedure.i QMulReal_and_Quaternion(R.f, *Q.Quaternion) 

 Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion))
 If *B
   *B\a=*Q\a * R: *B\b=*Q\b * R
   *B\c=*Q\c * R: *B\d=*Q\d * R
 EndIf
 ProcedureReturn *B

EndProcedure

Procedure.i QMulQuaternion(*Q1.Quaternion, *Q2.Quaternion) 

 Protected *B.Quaternion2=AllocateMemory(SizeOf(Quaternion2))
 If *B
   With *B
     \Qa\a= *Q1\a * *Q2\a: \Qa\b=-*Q1\b * *Q2\b
     \Qa\c=-*Q1\c * *Q2\c: \Qa\d=-*Q1\d * *Q2\d
     ;
     \Qb\a= *Q1\a * *Q2\b: \Qb\b= *Q1\b * *Q2\a
     \Qb\c= *Q1\c * *Q2\d: \Qb\d=-*Q1\d * *Q2\c
     ;
     \Qc\a= *Q1\a * *Q2\c: \Qc\b=-*Q1\b * *Q2\d
     \Qc\c= *Q1\c * *Q2\a: \Qc\d= *Q1\d * *Q2\b
     ;
     \Qd\a= *Q1\a * *Q2\d: \Qd\b= *Q1\b * *Q2\c
     \Qd\c=-*Q1\c * *Q2\b: \Qd\d= *Q1\d * *Q2\a
   EndWith
 EndIf
 ProcedureReturn *B

EndProcedure</lang> Implementation & test <lang PureBasic>Macro Show(QQ,NN=0) 

 StrF(QQ\a,NN)+","+StrF(QQ\b,NN)+","+StrF(QQ\c,NN)+","+StrF(QQ\d,NN)

EndMacro

Define.Quaternion A, B, C Define.f r=7 Define *X.Quaternion, *Y.Quaternion2

A\a=1: A\b=2: A\c=3: A\d=4 B\a=2: B\b=3: B\c=4: B\d=5 C\a=3: C\b=4: C\c=5: C\d=6

Debug "Q1= {"+Show(A,0)+"}" Debug "Q2= {"+Show(B,0)+"}" Debug "Q3= {"+Show(C,0)+"}"


Debug "Normal of Q1= "+StrF(QNorm(@A))

  • X=QNeg(@A)

Debug "Neg(Q1) ={"+Show(*X)+"}" : FreeMemory(*X)

  • X=QConj(@A)

Debug "Conj(Q1) ={"+Show(*X)+"}" : FreeMemory(*X)

  • X=QAddReal(r,@A)

Debug "r+A ={"+Show(*X)+"}" : FreeMemory(*X)

  • X=QAddQuaternion(@A,@B)

Debug "Q1+Q2 ={"+Show(*X)+"}" : FreeMemory(*X)

  • X=QAddQuaternion(@B,@C)

Debug "Q2+Q3 ={"+Show(*X)+"}" : FreeMemory(*X)

  • Y=QMulQuaternion(@A,@B)

Debug "Q1*Q2 =" Debug "{{"+Show(*Y\Qa)+"}" Debug " {"+Show(*Y\Qb)+"}" Debug " {"+Show(*Y\Qc)+"}" Debug " {"+Show(*Y\Qd)+"}}" FreeMemory(*Y)

  • Y=QMulQuaternion(@B,@A)

Debug "Q2*Q1 =" Debug "{{"+Show(*Y\Qa)+"}" Debug " {"+Show(*Y\Qb)+"}" Debug " {"+Show(*Y\Qc)+"}" Debug " {"+Show(*Y\Qd)+"}}" FreeMemory(*Y)</lang> Result 

Q1= {1,2,3,4}
Q2= {2,3,4,5}
Q3= {3,4,5,6}
Normal of Q1= 5.4772257805
Neg(Q1) ={-1,-2,-3,-4}
Conj(Q1) ={1,-2,-3,-4}
r+A      ={8,-2,-3,-4}
Q1+Q2    ={3,5,7,9}
Q2+Q3    ={5,7,9,11}
Q1*Q2    =
{{2,-6,-12,-20}
 {3,4,15,-16}
 {4,-10,6,12}
 {5,8,-9,8}}
Q2*Q1    =
{{2,-6,-12,-20}
 {4,3,16,-15}
 {6,-12,4,10}
 {8,9,-8,5}}

Python

This example extends Pythons namedtuples to add extra functionality. <lang python>from collections import namedtuple import math

class Q(namedtuple('Quaternion', 'real, i, j, k')): 

   'Quaternion type: Q(real=0.0, i=0.0, j=0.0, k=0.0)' 

   __slots__ = () 

   def __new__(_cls, real=0.0, i=0.0, j=0.0, k=0.0):
       'Defaults all parts of quaternion to zero'
       return super().__new__(_cls, float(real), float(i), float(j), float(k))

   def conjugate(self):
       return Q(self.real, -self.i, -self.j, -self.k)

   def _norm2(self):
       return sum( x*x for x in self)

   def norm(self):
       return math.sqrt(self._norm2())

   def reciprocal(self):
       n2 = self._norm2()
       return Q(*(x / n2 for x in self.conjugate())) 

   def __str__(self):
       'Shorter form of Quaternion as string'
       return 'Q(%g, %g, %g, %g)' % self

   def __neg__(self):
       return Q(-self.real, -self.i, -self.j, -self.k)

   def __add__(self, other):
       if type(other) == Q:
           return Q( *(s+o for s,o in zip(self, other)) )
       try:
           f = float(other)
       except:
           return NotImplemented
       return Q(self.real + f, self.i, self.j, self.k)

   def __radd__(self, other):
       return Q.__add__(self, other)

   def __mul__(self, other):
       if type(other) == Q:
           a1,b1,c1,d1 = self
           a2,b2,c2,d2 = other
           return Q(
                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 )
       try:
           f = float(other)
       except:
           return NotImplemented
       return Q(self.real * f, self.i * f, self.j * f, self.k * f)

   def __rmul__(self, other):
       return Q.__mul__(self, other)

   def __truediv__(self, other):
       if type(other) == Q:
           return self.__mul__(other.reciprocal())
       try:
           f = float(other)
       except:
           return NotImplemented
       return Q(self.real / f, self.i / f, self.j / f, self.k / f)

   def __rtruediv__(self, other):
       return other * self.reciprocal()

   __div__, __rdiv__ = __truediv__, __rtruediv__

Quaternion = Q

q = Q(1, 2, 3, 4) q1 = Q(2, 3, 4, 5) q2 = Q(3, 4, 5, 6) r = 7</lang>

Continued shell session Run the above with the -i flag to python on the command line, or run with idle then continue in the shell as follows: <lang python>>>> q Quaternion(real=1.0, i=2.0, j=3.0, k=4.0) >>> q1 Quaternion(real=2.0, i=3.0, j=4.0, k=5.0) >>> q2 Quaternion(real=3.0, i=4.0, j=5.0, k=6.0) >>> r 7 >>> q.norm() 5.477225575051661 >>> q1.norm() 7.3484692283495345 >>> q2.norm() 9.273618495495704 >>> -q Quaternion(real=-1.0, i=-2.0, j=-3.0, k=-4.0) >>> q.conjugate() Quaternion(real=1.0, i=-2.0, j=-3.0, k=-4.0) >>> r + q Quaternion(real=8.0, i=2.0, j=3.0, k=4.0) >>> q + r Quaternion(real=8.0, i=2.0, j=3.0, k=4.0) >>> q1 + q2 Quaternion(real=5.0, i=7.0, j=9.0, k=11.0) >>> q2 + q1 Quaternion(real=5.0, i=7.0, j=9.0, k=11.0) >>> q * r Quaternion(real=7.0, i=14.0, j=21.0, k=28.0) >>> r * q Quaternion(real=7.0, i=14.0, j=21.0, k=28.0) >>> q1 * q2 Quaternion(real=-56.0, i=16.0, j=24.0, k=26.0) >>> q2 * q1 Quaternion(real=-56.0, i=18.0, j=20.0, k=28.0) >>> assert q1 * q2 != q2 * q1 >>> >>> i, j, k = Q(0,1,0,0), Q(0,0,1,0), Q(0,0,0,1) >>> i*i Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> j*j Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> k*k Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> i*j*k Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> q1 / q2 Quaternion(real=0.7906976744186047, i=0.023255813953488358, j=-2.7755575615628914e-17, k=0.046511627906976744) >>> q1 / q2 * q2 Quaternion(real=2.0000000000000004, i=3.0000000000000004, j=4.000000000000001, k=5.000000000000001) >>> q2 * q1 / q2 Quaternion(real=2.0, i=3.465116279069768, j=3.906976744186047, k=4.767441860465116) >>> q1.reciprocal() * q1 Quaternion(real=0.9999999999999999, i=0.0, j=0.0, k=0.0) >>> q1 * q1.reciprocal() Quaternion(real=0.9999999999999999, i=0.0, j=0.0, k=0.0) >>> </lang>

Ruby

<lang ruby>class Quaternion 

 attr_reader :a, :b, :c, :d

 def initialize a, b, c, d
   @a, @b, @c, @d = a, b, c, d
 end

 def norm
   Math.sqrt(@a * @a + @b * @b + @c * @c + @d * @d)
 end

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

 def conj
   Quaternion.new(@a, -@b, -@c, -@d)
 end

 def + q
   if q.is_a? Quaternion
     Quaternion.new(@a + q.a, @b + q.b, @c + q.c, @d + q.d)
   elsif q.is_a? Numeric
     Quaternion.new(@a + q, @b, @c, @d)
   else
     q + self
   end
 end

 def - q
   self + -q
 end

 def * q
   if q.is_a? Quaternion
     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)
   elsif q.is_a? Numeric
     Quaternion.new(@a * q, @b * q, @c * q, @d * q)
   else
     q * self
   end
 end

 def == q
   @a == q.a and @b == q.b and @c == q.c and @d == q.d
 end

 def to_s
   "(#{@a}, #{@b}, #{@c}, #{@d})"
 end

 # this is Ruby's way of handling commutativity
 def coerce other
   return self, other
 end

end</lang>

IRB session 

>> require 'quaternions'
=> true
>> q = Quaternion.new(1, 2, 3, 4)
=> #<Quaternion:0x7fcd33128a88 @b=2, @a=1, @d=4, @c=3>
>> q1 = Quaternion.new(2, 3, 4, 5)
=> #<Quaternion:0x7fcd33122e80 @b=3, @a=2, @d=5, @c=4>
>> q2 = Quaternion.new(3, 4, 5, 6)
=> #<Quaternion:0x7fcd3311d278 @b=4, @a=3, @d=6, @c=5>
>> r = 7
=> 7
>> q.norm
=> 5.47722557505166
>> q1.norm
=> 7.34846922834953
>> q2.norm
=> 9.2736184954957
>> -q
=> #<Quaternion:0x7fcd331149c0 @b=-2, @a=-1, @d=-4, @c=-3>
>> q.conj
=> #<Quaternion:0x7fcd33112e18 @b=-2, @a=1, @d=-4, @c=-3>
>> q1 + q2
=> #<Quaternion:0x7fcd33110af0 @b=7, @a=5, @d=11, @c=9>
>> q2 + q1
=> #<Quaternion:0x7fcd3310e7c8 @b=7, @a=5, @d=11, @c=9>
>> q + r
=> #<Quaternion:0x7fcd3310c590 @b=2, @a=8, @d=4, @c=3>
>> r + q
=> #<Quaternion:0x7fcd3310a330 @b=2, @a=8, @d=4, @c=3>
>> q * r
=> #<Quaternion:0x7fcd33108170 @b=14, @a=7, @d=28, @c=21>
>> r * q
=> #<Quaternion:0x7fcd33105f88 @b=14, @a=7, @d=28, @c=21>
>> q1 * q2
=> #<Quaternion:0x7fcd33103cd8 @b=-50, @a=-56, @d=-34, @c=-40>
>> q2 * q1
=> #<Quaternion:0x7fcd33101a28 @b=-48, @a=-56, @d=-28, @c=-36>
>> q1 * q2 != q2 * q1
=> true

Tcl

Works with: Tcl version 8.6

or

Library: TclOO

<lang tcl>package require TclOO

  1. Support class that provides C++-like RAII lifetimes

oo::class create RAII-support { 

   constructor {} {

upvar 1 { end } end lappend end [self] trace add variable end unset [namespace code {my destroy}] 

   }
   destructor {

catch { upvar 1 { end } end trace remove variable end unset [namespace code {my destroy}] } 

   }
   method return Template:Level 1 {

incr level upvar 1 { end } end upvar $level { end } parent trace remove variable end unset [namespace code {my destroy}] lappend parent [self] trace add variable parent unset [namespace code {my destroy}] return -level $level [self] 

   }

}

  1. Class of quaternions

oo::class create Q { 

   superclass RAII-support
   variable R I J K
   constructor {{real 0} {i 0} {j 0} {k 0}} {

next namespace import ::tcl::mathfunc::* ::tcl::mathop::* variable R [double $real] I [double $i] J [double $j] K [double $k] 

   }
   self method return args {

[my new {*}$args] return 2 

   }

   method p {} {

return "Q($R,$I,$J,$K)" 

   }
   method values {} {

list $R $I $J $K 

   }

   method Norm {} {

+ [* $R $R] [* $I $I] [* $J $J] [* $K $K] 

   }

   method conjugate {} {

Q return $R [- $I] [- $J] [- $K] 

   }
   method norm {} {

sqrt [my Norm] 

   }
   method unit {} {

set n [my norm] Q return [/ $R $n] [/ $I $n] [/ $J $n] [/ $K $n] 

   }
   method reciprocal {} {

set n2 [my Norm] Q return [/ $R $n2] [/ $I $n2] [/ $J $n2] [/ $K $n2] 

   }
   method - Template:Q "" {

if {[llength [info level 0]] == 2} { Q return [- $R] [- $I] [- $J] [- $K] } [my + [$q -]] return 

   }
   method + q {

if {[info object isa object $q]} { lassign [$q values] real i j k Q return [+ $R $real] [+ $I $i] [+ $J $j] [+ $K $k] } Q return [+ $R [double $q]] $I $J $K 

   }
   method * q {

if {[info object isa object $q]} { lassign [my values] a1 b1 c1 d1 lassign [$q values] a2 b2 c2 d2 Q return [expr {$a1*$a2 - $b1*$b2 - $c1*$c2 - $d1*$d2}] \ [expr {$a1*$b2 + $b1*$a2 + $c1*$d2 - $d1*$c2}] \ [expr {$a1*$c2 - $b1*$d2 + $c1*$a2 + $d1*$b2}] \ [expr {$a1*$d2 + $b1*$c2 - $c1*$b2 + $d1*$a2}] } set f [double $q] Q return [* $R $f] [* $I $f] [* $J $f] [* $K $f] 

   }
   method == q {

expr { [info object isa object $q] && [info object isa typeof $q [self class]] && [my values] eq [$q values] } 

   }

   export - + * ==

}</lang> Demonstration code: <lang tcl>set q [Q new 1 2 3 4] set q1 [Q new 2 3 4 5] set q2 [Q new 3 4 5 6] set r 7

puts "q = [$q p]" puts "q1 = [$q1 p]" puts "q2 = [$q2 p]" puts "r = $r" puts "q norm = [$q norm]" puts "q1 norm = [$q1 norm]" puts "q2 norm = [$q2 norm]" puts "-q = [[$q -] p]" puts "q conj = [[$q conjugate] p]" puts "q + r = [[$q + $r] p]"

  1. Real numbers are not objects, so no extending operations for them

puts "q1 + q2 = [[$q1 + $q2] p]" puts "q2 + q1 = [[$q2 + $q1] p]" puts "q * r = [[$q * $r] p]" puts "q1 * q2 = [[$q1 * $q2] p]" puts "q2 * q1 = [[$q2 * $q1] p]" puts "equal(q1*q2, q2*q1) = [[$q1 * $q2] == [$q2 * $q1]]"</lang> 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
q norm = 5.477225575051661
q1 norm = 7.3484692283495345
q2 norm = 9.273618495495704
-q = Q(-1.0,-2.0,-3.0,-4.0)
q conj = 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)
equal(q1*q2, q2*q1) = 0
Retrieved from "https://rosettacode.org/wiki/Quaternion_type?oldid=153954"
Cookies help us deliver our services. By using our services, you agree to our use of cookies.
More information