Quaternion type

From Rosetta Code
Revision as of 18:19, 3 September 2010 by rosettacode>Blue Prawn
Quaternion type is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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

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>

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>

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>

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=153930"