# Quaternion type

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)

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:
${\displaystyle ={\sqrt {a^{2}+b^{2}+c^{2}+d^{2}}}}$
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)
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.

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

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:
# 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_print(r);
printf("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

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

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

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)

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

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

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

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