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:
- The norm of a quaternion:
- The negative of a quaternion:
=(-a, -b, -c, -d) - The conjugate of a quaternion:
=( a, -b, -c, -d) - Addition of a real number r and a quaternion q:
r + q = q + r = (a+r, b, c, d) - Addition of two quaternions:
q1 + q2 = (a1+a2, b1+b2, c1+c2, d1+d2) - Multiplication of a real number and a quaternion:
qr = rq = (ar, br, cr, dr) - 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 ) - Show that, for the two quaternions q1 and q2:
q1q2 != q2q1
If your language has built-in support for quaternions then use it.
C.f.
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
- note: This specimen retains the original python coding style.
<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);
- 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);
- 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);
- 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);
- 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);
- 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);
- 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>
- include <stdlib.h>
- include <stdbool.h>
- include <math.h>
typedef struct quaternion {
double q[4];
} quaternion_t;
quaternion_t *quaternion_new(void)
{
return malloc(sizeof(quaternion_t));
}
quaternion_t *quaternion_new_set(double q1, double q2, double q3, double q4) {
quaternion_t *q = malloc(sizeof(quaternion_t));
if (q != NULL) {
q->q[0] = q1; q->q[1] = q2; q->q[2] = q3; q->q[3] = q4;
}
return q;
}
void quaternion_copy(quaternion_t *r, quaternion_t *q)
{
size_t i;
if (r == NULL || q == NULL) return; for(i = 0; i < 4; i++) r->q[i] = q->q[i];
}
double quaternion_norm(quaternion_t *q)
{
size_t i;
double r = 0.0;
if (q == NULL) {
fprintf(stderr, "NULL quaternion in norm\n");
return 0.0;
}
for(i = 0; i < 4; i++) r += q->q[i] * q->q[i]; return sqrt(r);
}
void quaternion_neg(quaternion_t *r, quaternion_t *q)
{
size_t i;
if (q == NULL || r == NULL) return; for(i = 0; i < 4; i++) r->q[i] = -q->q[i];
}
void quaternion_conj(quaternion_t *r, quaternion_t *q)
{
size_t i;
if (q == NULL || r == NULL) return; r->q[0] = q->q[0]; for(i = 1; i < 4; i++) r->q[i] = -q->q[i];
}
void quaternion_add_d(quaternion_t *r, quaternion_t *q, double d)
{
if (q == NULL || r == NULL) return; quaternion_copy(r, q); r->q[0] += d;
}
void quaternion_add(quaternion_t *r, quaternion_t *a, quaternion_t *b)
{
size_t i;
if (r == NULL || a == NULL || b == NULL) return; for(i = 0; i < 4; i++) r->q[i] = a->q[i] + b->q[i];
}
void quaternion_mul_d(quaternion_t *r, quaternion_t *q, double d)
{
size_t i;
if (r == NULL || q == NULL) return; for(i = 0; i < 4; i++) r->q[i] = q->q[i] * d;
}
bool quaternion_equal(quaternion_t *a, quaternion_t *b) {
size_t i; for(i = 0; i < 4; i++) if (a->q[i] != b->q[i]) return false; return true;
}
- define A(N) (a->q[(N)])
- define B(N) (b->q[(N)])
- define R(N) (r->q[(N)])
void quaternion_mul(quaternion_t *r, quaternion_t *a, quaternion_t *b) {
size_t i; double ri = 0.0;
if (r == NULL || a == NULL || b == NULL) return; R(0) = A(0)*B(0) - A(1)*B(1) - A(2)*B(2) - A(3)*B(3); R(1) = A(0)*B(1) + A(1)*B(0) + A(2)*B(3) - A(3)*B(2); R(2) = A(0)*B(2) - A(1)*B(3) + A(2)*B(0) + A(3)*B(1); R(3) = A(0)*B(3) + A(1)*B(2) - A(2)*B(1) + A(3)*B(0);
}
- undef A
- undef B
- undef R
void quaternion_print(quaternion_t *q)
{
if (q == NULL) return;
printf("(%lf, %lf, %lf, %lf)\n",
q->q[0], q->q[1], q->q[2], q->q[3]); }</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++
This example uses templates to provide the underlying data-type, and includes several extra functions and constructors that often come up when using quaternions.
<lang cpp>#include <iostream> using namespace std;
template<class T = double> class Quaternion { public:
T w, x, y, z;
// Numerical constructor
Quaternion(const T &w, const T &x, const T &y, const T &z): w(w), x(x), y(y), z(z) {};
Quaternion(const T &x, const T &y, const T &z): w(T()), x(x), y(y), z(z) {}; // For 3-rotations
Quaternion(const T &r): w(r), x(T()), y(T()), z(T()) {};
Quaternion(): w(T()), x(T()), y(T()), z(T()) {};
// Copy constructor and assignment
Quaternion(const Quaternion &q): w(q.w), x(q.x), y(q.y), z(q.z) {};
Quaternion& operator=(const Quaternion &q) { w=q.w; x=q.x; y=q.y; z=q.z; return *this; }
// Unary operators
Quaternion operator-() const { return Quaternion(-w, -x, -y, -z); }
Quaternion operator~() const { return Quaternion(w, -x, -y, -z); } // Conjugate
// Norm-squared. SQRT would have to be made generic to be used here
T normSquared() const { return w*w + x*x + y*y + z*z; }
// In-place operators
Quaternion& operator+=(const T &r)
{ w += r; return *this; }
Quaternion& operator+=(const Quaternion &q)
{ w += q.w; x += q.x; y += q.y; z += q.z; return *this; }
Quaternion& operator-=(const T &r)
{ w -= r; return *this; }
Quaternion& operator-=(const Quaternion &q)
{ w -= q.w; x -= q.x; y -= q.y; z -= q.z; return *this; }
Quaternion& operator*=(const T &r)
{ w *= r; x *= r; y *= r; z *= r; return *this; }
Quaternion& operator*=(const Quaternion &q)
{
T oldW(w), oldX(x), oldY(y), oldZ(z);
w = oldW*q.w - oldX*q.x - oldY*q.y - oldZ*q.z;
x = oldW*q.x + oldX*q.w + oldY*q.z - oldZ*q.y;
y = oldW*q.y + oldY*q.w + oldZ*q.x - oldX*q.z;
z = oldW*q.z + oldZ*q.w + oldX*q.y - oldY*q.x;
return *this;
}
Quaternion& operator/=(const T &r)
{ w /= r; x /= r; y /= r; z /= r; return *this; }
Quaternion& operator/=(const Quaternion &q)
{
T oldW(w), oldX(x), oldY(y), oldZ(z), n(q.normSquared());
w = (oldW*q.w + oldX*q.x + oldY*q.y + oldZ*q.z) / n;
x = (oldX*q.w - oldW*q.x + oldY*q.z - oldZ*q.y) / n;
y = (oldY*q.w - oldW*q.y + oldZ*q.x - oldX*q.z) / n;
z = (oldZ*q.w - oldW*q.z + oldX*q.y - oldY*q.x) / n;
return *this;
}
// Binary operators based on in-place operators
Quaternion operator+(const T &r) const { return Quaternion(*this) += r; }
Quaternion operator+(const Quaternion &q) const { return Quaternion(*this) += q; }
Quaternion operator-(const T &r) const { return Quaternion(*this) -= r; }
Quaternion operator-(const Quaternion &q) const { return Quaternion(*this) -= q; }
Quaternion operator*(const T &r) const { return Quaternion(*this) *= r; }
Quaternion operator*(const Quaternion &q) const { return Quaternion(*this) *= q; }
Quaternion operator/(const T &r) const { return Quaternion(*this) /= r; }
Quaternion operator/(const Quaternion &q) const { return Quaternion(*this) /= q; }
// Comparison operators, as much as they make sense
bool operator==(const Quaternion &q) const
{ return (w == q.w) && (x == q.x) && (y == q.y) && (z == q.z); }
bool operator!=(const Quaternion &q) const { return !operator==(q); }
// The operators above allow quaternion op real. These allow real op quaternion. // Uses the above where appropriate. template<class T> friend Quaternion<T> operator+(const T &r, const Quaternion<T> &q); template<class T> friend Quaternion<T> operator-(const T &r, const Quaternion<T> &q); template<class T> friend Quaternion<T> operator*(const T &r, const Quaternion<T> &q); template<class T> friend Quaternion<T> operator/(const T &r, const Quaternion<T> &q); // Allows cout << q template<class T> friend ostream& operator<<(ostream &io, const Quaternion<T> &q);
};
// Friend functions need to be outside the actual class definition template<class T> Quaternion<T> operator+(const T &r, const Quaternion<T> &q)
{ return q+r; }
template<class T> Quaternion<T> operator-(const T &r, const Quaternion<T> &q)
{ return Quaternion<T>(r-q.w, q.x, q.y, q.z); }
template<class T> Quaternion<T> operator*(const T &r, const Quaternion<T> &q)
{ return q*r; }
template<class T> Quaternion<T> operator/(const T &r, const Quaternion<T> &q) {
T n(q.normSquared()); return Quaternion(r*q.w/n, -r*q.x/n, -r*q.y/n, -r*q.z/n);
}
template<class T> ostream& operator<<(ostream &io, const Quaternion<T> &q) {
io << q.w; (q.x < T()) ? (io << " - " << (-q.x) << "i") : (io << " + " << q.x << "i"); (q.y < T()) ? (io << " - " << (-q.y) << "j") : (io << " + " << q.y << "j"); (q.z < T()) ? (io << " - " << (-q.z) << "k") : (io << " + " << q.z << "k"); return io;
}</lang>
Test program: <lang cpp>int main() {
Quaternion<> q0(1, 2, 3, 4); Quaternion<> q1(2, 3, 4, 5); Quaternion<> q2(3, 4, 5, 6); double r = 7;
cout << "q0: " << q0 << endl; cout << "q1: " << q1 << endl; cout << "q2: " << q2 << endl; cout << "r: " << r << endl; cout << endl; cout << "-q0: " << -q0 << endl; cout << "~q0: " << ~q0 << endl; cout << endl; cout << "r * q0: " << r*q0 << endl; cout << "r + q0: " << r+q0 << endl; cout << "q0 / r: " << q0/r << endl; cout << "q0 - r: " << q0-r << endl; cout << endl; cout << "q0 + q1: " << q0+q1 << endl; cout << "q0 - q1: " << q0-q1 << endl; cout << "q0 * q1: " << q0*q1 << endl; cout << "q0 / q1: " << q0/q1 << endl; cout << endl; cout << "q0 * ~q0: " << q0*~q0 << endl; cout << "q0 + q1*q2: " << q0+q1*q2 << endl; cout << "(q0 + q1)*q2: " << (q0+q1)*q2 << endl; cout << "q0*q1*q2: " << q0*q1*q2 << endl; cout << "(q0*q1)*q2: " << (q0*q1)*q2 << endl; cout << "q0*(q1*q2): " << q0*(q1*q2) << endl; cout << endl; cout << "||q0||: " << sqrt(q0.normSquared()) << endl; cout << endl; cout << "q0*q1 - q1*q0: " << (q0*q1 - q1*q0) << endl;
// Other base types Quaternion<int> q5(2), q6(3); cout << endl << q5*q6 << endl;
}</lang>
Output:
q0: 1 + 2i + 3j + 4k q1: 2 + 3i + 4j + 5k q2: 3 + 4i + 5j + 6k r: 7 -q0: -1 - 2i - 3j - 4k ~q0: 1 - 2i - 3j - 4k r * q0: 7 + 14i + 21j + 28k r + q0: 8 + 2i + 3j + 4k q0 / r: 0.142857 + 0.285714i + 0.428571j + 0.571429k q0 - r: -6 + 2i + 3j + 4k q0 + q1: 3 + 5i + 7j + 9k q0 - q1: -1 - 1i - 1j - 1k q0 * q1: -36 + 6i + 12j + 12k q0 / q1: 0.740741 + 0i + 0.0740741j + 0.037037k q0 * ~q0: 30 + 0i + 0j + 0k q0 + q1*q2: -55 + 18i + 27j + 30k (q0 + q1)*q2: -100 + 24i + 42j + 42k q0*q1*q2: -264 - 114i - 132j - 198k (q0*q1)*q2: -264 - 114i - 132j - 198k q0*(q1*q2): -264 - 114i - 132j - 198k ||q0||: 5.47723 q0*q1 - q1*q0: 0 - 2i + 4j - 2k 6 + 0i + 0j + 0k
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
<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>
E
<lang e>interface Quaternion guards QS {} def makeQuaternion(a, b, c, d) {
return def quaternion implements QS {
to __printOn(out) {
out.print("(", a, " + ", b, "i + ")
out.print(c, "j + ", d, "k)")
}
# Task requirement 1
to norm() {
return (a**2 + b**2 + c**2 + d**2).sqrt()
}
# Task requirement 2
to negate() {
return makeQuaternion(-a, -b, -c, -d)
}
# Task requirement 3
to conjugate() {
return makeQuaternion(a, -b, -c, -d)
}
# Task requirement 4, 5
# This implements q + r; r + q is deliberately prohibited by E
to add(other :any[Quaternion, int, float64]) {
switch (other) {
match q :Quaternion {
return makeQuaternion(
a+q.a(), b+q.b(), c+q.c(), d+q.d())
}
match real {
return makeQuaternion(a+real, b, c, d)
}
}
}
# Task requirement 6, 7
# This implements q * r; r * q is deliberately prohibited by E
to multiply(other :any[Quaternion, int, float64]) {
switch (other) {
match q :Quaternion {
return makeQuaternion(
a*q.a() - b*q.b() - c*q.c() - d*q.d(),
a*q.b() + b*q.a() + c*q.d() - d*q.c(),
a*q.c() - b*q.d() + c*q.a() + d*q.b(),
a*q.d() + b*q.c() - c*q.b() + d*q.a())
}
match real {
return makeQuaternion(real*a, real*b, real*c, real*d)
}
}
}
to a() { return a }
to b() { return b }
to c() { return c }
to d() { return d }
}
}</lang>
<lang e>? def q1 := makeQuaternion(2,3,4,5)
- value: (2 + 3i + 4j + 5k)
? def q2 := makeQuaternion(3,4,5,6)
- value: (3 + 4i + 5j + 6k)
? q1+q2
- value: (5 + 7i + 9j + 11k)
? q1*q2
- value: (-56 + 16i + 24j + 26k)
? q2*q1
- value: (-56 + 18i + 20j + 28k)
? q1+(-2)
- value: (0 + 3i + 4j + 5k)</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
<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
GAP
<lang gap># GAP has built-in support for quaternions
A := QuaternionAlgebra(Rationals);
- <algebra-with-one of dimension 4 over Rationals>
b := BasisVectors(Basis(A));
- [ e, i, j, k ]
q := [1, 2, 3, 4]*b;
- e+(2)*i+(3)*j+(4)*k
- Computing norm may be difficult, since the result would be in a quadratic field (Sqrt exists in GAP, is quite unusual, see ?E in GAP, and the following example
Sqrt(5/3);
- 1/3*E(60)^7+1/3*E(60)^11-1/3*E(60)^19-1/3*E(60)^23-1/3*E(60)^31+1/3*E(60)^43-1/3*E(60)^47+1/3*E(60)^59
q1 := [2, 3, 4, 5]*b;
- (2)*e+(3)*i+(4)*j+(5)*k
q2 := [3, 4, 5, 6]*b;
- (3)*e+(4)*i+(5)*j+(6)*k
q1*q2 - q2*q1;
- (-2)*i+(4)*j+(-2)*k
- Can't add directly to a rational, one must make a quaternion of it
r:=5/3*b[1];
- (5/3)*e
r + q;
- (8/3)*e+(2)*i+(3)*j+(4)*k
- For multiplication, no problem (we are in an algebra over rationals !)
r * q;
- (5/3)*e+(10/3)*i+(5)*j+(20/3)*k
5/3 * q;
- (5/3)*e+(10/3)*i+(5)*j+(20/3)*k
- Negative
-q; (-1)*e+(-2)*i+(-3)*j+(-4)*k</lang>
Go
Conventions for method receiver, parameter, and return values modeled after Go's big number package. It provides flexibility without requiring unnecessary object creation. The test program creates only four quaternion objects, the three inputs and one more for an output. The three inputs are reused repeatedly without being modified. The output is also reused repeatedly, being overwritten for each operation. <lang go>package main
import (
"fmt" "math"
)
type qtn struct {
r, i, j, k float64
}
var (
q = &qtn{1, 2, 3, 4}
q1 = &qtn{2, 3, 4, 5}
q2 = &qtn{3, 4, 5, 6}
r float64 = 7
)
func main() {
fmt.Println("Inputs")
fmt.Println("q:", q)
fmt.Println("q1:", q1)
fmt.Println("q2:", q2)
fmt.Println("r:", r)
fmt.Println("\nFunctions")
qr := &qtn{}
fmt.Println("q.norm():", q.norm())
fmt.Println("neg(q):", qr.neg(q))
fmt.Println("conj(q):", qr.conj(q))
fmt.Println("addF(q, r):", qr.addF(q, r))
fmt.Println("addQ(q1, q2):", qr.addQ(q1, q2))
fmt.Println("mulF(q, r):", qr.mulF(q, r))
fmt.Println("mulQ(q1, q2):", qr.mulQ(q1, q2))
fmt.Println("mulQ(q2, q1):", qr.mulQ(q2, q1))
}
func (q *qtn) String() string {
return fmt.Sprintf("(%g, %g, %g, %g)", q.r, q.i, q.j, q.k)
}
func (q *qtn) norm() float64 {
return math.Sqrt(q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k)
}
func (z *qtn) neg(q *qtn) *qtn {
z.r, z.i, z.j, z.k = -q.r, -q.i, -q.j, -q.k return z
}
func (z *qtn) conj(q *qtn) *qtn {
z.r, z.i, z.j, z.k = q.r, -q.i, -q.j, -q.k return z
}
func (z *qtn) addF(q *qtn, r float64) *qtn {
z.r, z.i, z.j, z.k = q.r+r, q.i, q.j, q.k return z
}
func (z *qtn) addQ(q1, q2 *qtn) *qtn {
z.r, z.i, z.j, z.k = q1.r+q2.r, q1.i+q2.i, q1.j+q2.j, q1.k+q2.k return z
}
func (z *qtn) mulF(q *qtn, r float64) *qtn {
z.r, z.i, z.j, z.k = q.r*r, q.i*r, q.j*r, q.k*r return z
}
func (z *qtn) mulQ(q1, q2 *qtn) *qtn {
z.r = q1.r*q2.r - q1.i*q2.i - q1.j*q2.j - q1.k*q2.k z.i = q1.r*q2.i + q1.i*q2.r + q1.j*q2.k - q1.k*q2.j z.j = q1.r*q2.j - q1.i*q2.k + q1.j*q2.r + q1.k*q2.i z.k = q1.r*q2.k + q1.i*q2.j - q1.j*q2.i + q1.k*q2.r return z
}</lang> Output:
Inputs q: (1, 2, 3, 4) q1: (2, 3, 4, 5) q2: (3, 4, 5, 6) r: 7 Functions q.norm(): 5.477225575051661 neg(q): (-1, -2, -3, -4) conj(q): (1, -2, -3, -4) addF(q, r): (8, 2, 3, 4) addQ(q1, q2): (5, 7, 9, 11) mulF(q, r): (7, 14, 21, 28) mulQ(q1, q2): (-56, 16, 24, 26) mulQ(q2, q1): (-56, 18, 20, 28)
Haskell
<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]
Icon and Unicon
Using Unicon's class system.
<lang Unicon> class Quaternion(a, b, c, d)
method norm () return sqrt (a*a + b*b + c*c + d*d) end
method negative () return Quaternion(-a, -b, -c, -d) end
method conjugate () return Quaternion(a, -b, -c, -d) end
method add (n)
if type(n) == "Quaternion__state"
then return Quaternion(a+n.a, b+n.b, c+n.c, d+n.d)
else return Quaternion(a+n, b, c, d)
end
method multiply (n)
if type(n) == "Quaternion__state"
then return Quaternion(a*n.a - b*n.b - c*n.c - d*n.d,
a*n.b + b*n.a + c*n.d - d*n.c,
a*n.c - b*n.d + c*n.a + d*n.b,
a*n.d + b*n.c - c*n.b + d*n.a)
else return Quaternion(a*n, b*n, c*n, d*n)
end
method sign (n) return if n >= 0 then "+" else "-" end
method string ()
return ("" || a || sign(b) || abs(b) || "i" || sign(c) || abs(c) || "j" || sign(d) || abs(d) || "k");
end
initially(a, b, c, d) self.a := if /a then 0 else a self.b := if /b then 0 else b self.c := if /c then 0 else c self.d := if /d then 0 else d
end </lang>
To test the above:
<lang Unicon> procedure main ()
q := Quaternion (1,2,3,4) q1 := Quaternion (2,3,4,5) q2 := Quaternion (3,4,5,6) r := 7
write ("The norm of " || q.string() || " is " || q.norm ())
write ("The negative of " || q.string() || " is " || q.negative().string ())
write ("The conjugate of " || q.string() || " is " || q.conjugate().string ())
write ("Sum of " || q.string() || " and " || r || " is " || q.add(r).string ())
write ("Sum of " || q.string() || " and " || q1.string() || " is " || q.add(q1).string ())
write ("Product of " || q.string() || " and " || r || " is " || q.multiply(r).string ())
write ("Product of " || q.string() || " and " || q1.string() || " is " || q.multiply(q1).string ())
write ("q1*q2 = " || q1.multiply(q2).string ())
write ("q2*q1 = " || q2.multiply(q1).string ())
end </lang>
Output:
The norm of 1+2i+3j+4k is 5.477225575 The negative of 1+2i+3j+4k is -1-2i-3j-4k The conjugate of 1+2i+3j+4k is 1-2i-3j-4k Sum of 1+2i+3j+4k and 7 is 8+2i+3j+4k Sum of 1+2i+3j+4k and 2+3i+4j+5k is 3+5i+7j+9k Product of 1+2i+3j+4k and 7 is 7+14i+21j+28k Product of 1+2i+3j+4k and 2+3i+4j+5k is -36+6i+12j+12k q1*q2 = -56+16i+24j+26k q2*q1 = -56+18i+20j+28k
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>
Java
<lang java>public class Quaternion {
private final double a, b, c, d;
public Quaternion(double a, double b, double c, double d) {
this.a = a;
this.b = b;
this.c = c;
this.d = d;
}
public Quaternion(double r) {
this(r, 0.0, 0.0, 0.0);
}
public double norm() {
return Math.sqrt(a * a + b * b + c * c + d * d);
}
public Quaternion negative() {
return new Quaternion(-a, -b, -c, -d);
}
public Quaternion conjugate() {
return new Quaternion(a, -b, -c, -d);
}
public Quaternion add(double r) {
return new Quaternion(a + r, b, c, d);
}
public static Quaternion add(Quaternion q, double r) {
return q.add(r);
}
public static Quaternion add(double r, Quaternion q) {
return q.add(r);
}
public Quaternion add(Quaternion q) {
return new Quaternion(a + q.a, b + q.b, c + q.c, d + q.d);
}
public static Quaternion add(Quaternion q1, Quaternion q2) {
return q1.add(q2);
}
public Quaternion times(double r) {
return new Quaternion(a * r, b * r, c * r, d * r);
}
public static Quaternion times(Quaternion q, double r) {
return q.times(r);
}
public static Quaternion times(double r, Quaternion q) {
return q.times(r);
}
public Quaternion times(Quaternion q) {
return new Quaternion(
a * q.a - b * q.b - c * q.c - d * q.d,
a * q.b + b * q.a + c * q.d - d * q.c,
a * q.c - b * q.d + c * q.a + d * q.b,
a * q.d + b * q.c - c * q.b + d * q.a
);
}
public static Quaternion times(Quaternion q1, Quaternion q2) {
return q1.times(q2);
}
@Override
public boolean equals(Object obj) {
if (!(obj instanceof Quaternion)) return false;
final Quaternion other = (Quaternion) obj;
if (Double.doubleToLongBits(this.a) != Double.doubleToLongBits(other.a)) return false;
if (Double.doubleToLongBits(this.b) != Double.doubleToLongBits(other.b)) return false;
if (Double.doubleToLongBits(this.c) != Double.doubleToLongBits(other.c)) return false;
if (Double.doubleToLongBits(this.d) != Double.doubleToLongBits(other.d)) return false;
return true;
}
@Override
public String toString() {
return String.format("%.2f + %.2fi + %.2fj + %.2fk", a, b, c, d).replaceAll("\\+ -", "- ");
}
public String toQuadruple() {
return String.format("(%.2f, %.2f, %.2f, %.2f)", a, b, c, d);
}
public static void main(String[] args) {
Quaternion q = new Quaternion(1.0, 2.0, 3.0, 4.0);
Quaternion q1 = new Quaternion(2.0, 3.0, 4.0, 5.0);
Quaternion q2 = new Quaternion(3.0, 4.0, 5.0, 6.0);
double r = 7.0;
System.out.format("q = %s%n", q);
System.out.format("q1 = %s%n", q1);
System.out.format("q2 = %s%n", q2);
System.out.format("r = %.2f%n%n", r);
System.out.format("\u2016q\u2016 = %.2f%n", q.norm());
System.out.format("-q = %s%n", q.negative());
System.out.format("q* = %s%n", q.conjugate());
System.out.format("q + r = %s%n", q.add(r));
System.out.format("q1 + q2 = %s%n", q1.add(q2));
System.out.format("q \u00d7 r = %s%n", q.times(r));
Quaternion q1q2 = q1.times(q2);
Quaternion q2q1 = q2.times(q1);
System.out.format("q1 \u00d7 q2 = %s%n", q1q2);
System.out.format("q2 \u00d7 q1 = %s%n", q2q1);
System.out.format("q1 \u00d7 q2 %s q2 \u00d7 q1%n", (q1q2.equals(q2q1) ? "=" : "\u2260"));
}
}</lang>
This outputs:
q = 1.00 + 2.00i + 3.00j + 4.00k q1 = 2.00 + 3.00i + 4.00j + 5.00k q2 = 3.00 + 4.00i + 5.00j + 6.00k r = 7.00 ‖q‖ = 5.48 -q = -1.00 - 2.00i - 3.00j - 4.00k q* = 1.00 - 2.00i - 3.00j - 4.00k q + r = 8.00 + 2.00i + 3.00j + 4.00k q1 + q2 = 5.00 + 7.00i + 9.00j + 11.00k q × r = 7.00 + 14.00i + 21.00j + 28.00k q1 × q2 = -56.00 + 16.00i + 24.00j + 26.00k q2 × q1 = -56.00 + 18.00i + 20.00j + 28.00k q1 × q2 ≠ q2 × q1
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
PARI/GP
Here is a simple solution in GP. I think it's possible to implement this type directly in Pari by abusing t_COMPLEX, but I haven't attempted this. <lang>q.norm={ if(type(q) != "t_VEC" || #q != 4, error("incorrect type")); sqrt(q[1]^2+q[2]^2+q[3]^2+q[4]^2) }; q.conj={ if(type(q) != "t_VEC" || #q != 4, error("incorrect type")); -[-q[1],q[2],q[3],q[4]] }; q.add={ if(type(q) != "t_VEC" || #q != 4, error("incorrect type")); x->if(type(x) == "t_INT" || type(x) == t_REAL, [q[1]+x,q[2],q[3],q[4]] , if(type(x) == "t_VEC" && #x == 4, q+x , error("incorrect type") ) ) }; q.mult={ if(type(q) != "t_VEC" || #q != 4, error("incorrect type")); x->if(type(x) == "t_INT" || type(x) == t_REAL, x*q , if(type(x) == "t_VEC" && #x == 4, [q[1]*x[1] - q[2]*x[2] - q[3]*x[3] - q[4]*x[4], q[1]*x[2] + q[2]*x[1] + q[3]*x[4] - q[4]*x[3], q[1]*x[3] - q[2]*x[4] + q[3]*x[1] + q[4]*x[2], q[1]*x[4] + q[2]*x[3] - q[3]*x[2] + q[4]*x[1]] , error("incorrect type") ) ) };</lang> Usage: <lang>r=7;q=[1,2,3,4];q1=[2,3,4,5];q2=[3,4,5,6]; q.norm -q q.conj q.add(r) q1.add(q2) q1.add(q2) \\ or q1+q2 q.mult(r) \\ or r*q or q*r q1.mult(q2) q1.mult(q2) != q2.mult(q1)</lang>
Perl
<lang Perl>package Quaternion; use List::Util 'reduce'; use List::MoreUtils 'pairwise';
sub make {
my $cls = shift;
if (@_ == 1) { return bless [ @_, 0, 0, 0 ] }
elsif (@_ == 4) { return bless [ @_ ] }
else { die "Bad number of components: @_" }
}
sub _abs { sqrt reduce { $a + $b * $b } @{ +shift } } sub _neg { bless [ map(-$_, @{+shift}) ] } sub _str { "(@{+shift})" }
sub _add {
my ($x, $y) = @_;
$y = [ $y, 0, 0, 0 ] unless ref $y;
bless [ pairwise { $a + $b } @$x, @$y ]
}
sub _sub {
my ($x, $y, $swap) = @_;
$y = [ $y, 0, 0, 0 ] unless ref $y;
my @x = pairwise { $a - $b } @$x, @$y;
if ($swap) { $_ = -$_ for @x }
bless \@x;
}
sub _mul {
my ($x, $y) = @_;
if (!ref $y) { return bless [ map($_ * $y, @$x) ] }
my ($a1, $b1, $c1, $d1) = @$x;
my ($a2, $b2, $c2, $d2) = @$y;
bless [ $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]
}
sub conjugate {
my @a = map { -$_ } @{$_[0]};
$a[0] = $_[0][0];
bless \@a
}
use overload (
'""' => \&_str,
'+' => \&_add,
'-' => \&_sub,
'*' => \&_mul,
'neg' => \&_neg,
'abs' => \&_abs,
);
package main;
my $a = Quaternion->make(1, 2, 3, 4); my $b = Quaternion->make(1, 1, 1, 1);
print "a = $a\n"; print "b = $b\n"; print "|a| = ", abs($a), "\n"; print "-a = ", -$a, "\n"; print "a + 1 = ", $a + 1, "\n"; print "a + b = ", $a + $b, "\n"; print "a - b = ", $a - $b, "\n"; print "a conjugate is ", $a->conjugate, "\n"; print "a * b = ", $a * $b, "\n"; print "b * a = ", $b * $a, "\n";</lang>
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
PicoLisp
<lang PicoLisp>(scl 6)
(def 'quatCopy copy)
(de quatNorm (Q)
(sqrt (apply + (mapcar * Q Q))) )
(de quatNeg (Q)
(mapcar - Q) )
(de quatConj (Q)
(cons (car Q) (mapcar - (cdr Q))) )
(de quatAddR (Q R)
(cons (+ R (car Q)) (cdr Q)) )
(de quatAdd (Q1 Q2)
(mapcar + Q1 Q2) )
(de quatMulR (Q R)
(mapcar */ (mapcar * Q (circ R)) (1.0 .)) )
(de quatMul (Q1 Q2)
(mapcar
'((Ops I)
(sum '((Op R I) (Op (*/ R (get Q2 I) 1.0))) Ops Q1 I) )
'((+ - - -) (+ + + -) (+ - + +) (+ + - +))
'((1 2 3 4) (2 1 4 3) (3 4 1 2) (4 3 2 1)) ) )
(de quatFmt (Q)
(mapcar '((R S) (pack (format R *Scl) S))
Q
'(" + " "i + " "j + " "k") ) )</lang>
Test: <lang PicoLisp>(setq
Q (1.0 2.0 3.0 4.0) Q1 (2.0 3.0 4.0 5.0) Q2 (3.0 4.0 5.0 6.0) R 7.0 )
(prinl "R = " (format R *Scl)) (prinl "Q = " (quatFmt Q)) (prinl "Q1 = " (quatFmt Q1)) (prinl "Q2 = " (quatFmt Q2)) (prinl) (prinl "norm(Q) = " (format (quatNorm Q) *Scl)) (prinl "norm(Q1) = " (format (quatNorm Q1) *Scl)) (prinl "norm(Q2) = " (format (quatNorm Q2) *Scl)) (prinl "neg(Q) = " (quatFmt (quatNeg Q))) (prinl "conj(Q) = " (quatFmt (quatConj Q))) (prinl "Q + R = " (quatFmt (quatAddR Q R))) (prinl "Q1 + Q2 = " (quatFmt (quatAdd Q1 Q2))) (prinl "Q * R = " (quatFmt (quatMulR Q R))) (prinl "Q1 * Q2 = " (quatFmt (quatMul Q1 Q2))) (prinl "Q2 * Q1 = " (quatFmt (quatMul Q2 Q1))) (prinl (if (= (quatMul Q1 Q2) (quatMul Q2 Q1)) "Equal" "Not equal"))</lang> Output:
R = 7.000000 Q = 1.000000 + 2.000000i + 3.000000j + 4.000000k Q1 = 2.000000 + 3.000000i + 4.000000j + 5.000000k Q2 = 3.000000 + 4.000000i + 5.000000j + 6.000000k norm(Q) = 5.477225 norm(Q1) = 7.348469 norm(Q2) = 9.273618 neg(Q) = -1.000000 + -2.000000i + -3.000000j + -4.000000k conj(Q) = 1.000000 + -2.000000i + -3.000000j + -4.000000k Q + R = 8.000000 + 2.000000i + 3.000000j + 4.000000k Q1 + Q2 = 5.000000 + 7.000000i + 9.000000j + 11.000000k Q * R = 7.000000 + 14.000000i + 21.000000j + 28.000000k Q1 * Q2 = -56.000000 + 16.000000i + 24.000000j + 26.000000k Q2 * Q1 = -56.000000 + 18.000000i + 20.000000j + 28.000000k Not equal
PureBasic
<lang PureBasic>Structure Quaternion
a.f b.f c.f d.f
EndStructure
Procedure.f QNorm(*x.Quaternion)
ProcedureReturn Sqr(Pow(*x\a, 2) + Pow(*x\b, 2) + Pow(*x\c, 2) + Pow(*x\d, 2))
EndProcedure
- If supplied, the result is returned in the quaternion structure *res,
- otherwise a new quaternion is created. A pointer to the result is returned.
Procedure QNeg(*x.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf If *res *res\a = -*x\a *res\b = -*x\b *res\c = -*x\c *res\d = -*x\d EndIf ProcedureReturn *res
EndProcedure
Procedure QConj(*x.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf If *res *res\a = *x\a *res\b = -*x\b *res\c = -*x\c *res\d = -*x\d EndIf ProcedureReturn *res
EndProcedure
Procedure QAddReal(r.f, *x.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf If *res *res\a = *x\a + r *res\b = *x\b *res\c = *x\c *res\d = *x\d EndIf ProcedureReturn *res
EndProcedure
Procedure QAddQuaternion(*x.Quaternion, *y.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf If *res *res\a = *x\a + *y\a *res\b = *x\b + *y\b *res\c = *x\c + *y\c *res\d = *x\d + *y\d EndIf ProcedureReturn *res
EndProcedure
Procedure QMulReal_and_Quaternion(r.f, *x.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf If *res *res\a = *x\a * r *res\b = *x\b * r *res\c = *x\c * r *res\d = *x\d * r EndIf ProcedureReturn *res
EndProcedure
Procedure QMulQuaternion(*x.Quaternion, *y.Quaternion, *res.Quaternion = 0)
If *res = 0: *res.Quaternion = AllocateMemory(SizeOf(Quaternion)): EndIf
If *res
*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
EndIf
ProcedureReturn *res
EndProcedure
Procedure Q_areEqual(*x.Quaternion, *y.Quaternion)
If (*x\a <> *y\a) Or (*x\b <> *y\b) Or (*x\c <> *y\c) Or (*x\d <> *y\d) ProcedureReturn 0 ;false EndIf ProcedureReturn 1 ;true
EndProcedure</lang> Implementation & test <lang PureBasic>Procedure.s ShowQ(*x.Quaternion, NN = 0)
ProcedureReturn "{" + StrF(*x\a, NN) + "," + StrF(*x\b, NN) + "," + StrF(*x\c, NN) + "," + StrF(*x\d, NN) + "}"
EndProcedure
If OpenConsole()
Define.Quaternion Q0, Q1, Q2, res, res2
Define.f r = 7
Q0\a = 1: Q0\b = 2: Q0\c = 3: Q0\d = 4
Q1\a = 2: Q1\b = 3: Q1\c = 4: Q1\d = 5
Q2\a = 3: Q2\b = 4: Q2\c = 5: Q2\d = 6
PrintN("Q0 = " + ShowQ(Q0, 0))
PrintN("Q1 = " + ShowQ(Q1, 0))
PrintN("Q2 = " + ShowQ(Q2, 0))
PrintN("Normal of Q0 = " + StrF(QNorm(Q0)))
PrintN("Neg(Q0) = " + ShowQ(QNeg(Q0, res)))
PrintN("Conj(Q0) = " + ShowQ(QConj(Q0, res)))
PrintN("r + Q0 = " + ShowQ(QAddReal(r, Q0, res)))
PrintN("Q0 + Q1 = " + ShowQ(QAddQuaternion(Q0, Q1, res)))
PrintN("Q1 + Q2 = " + ShowQ(QAddQuaternion(Q1, Q2, res)))
PrintN("Q1 * Q2 = " + ShowQ(QMulQuaternion(Q1, Q2, res)))
PrintN("Q2 * Q1 = " + ShowQ(QMulQuaternion(Q2, Q1, res2)))
Print( "Q1 * Q2"): If Q_areEqual(res, res2): Print(" = "): Else: Print(" <> "): EndIf: Print( "Q2 * Q1")
Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
CloseConsole()
EndIf</lang> Result
Q0 = {1,2,3,4}
Q1 = {2,3,4,5}
Q2 = {3,4,5,6}
Normal of Q0 = 5.4772257805
Neg(Q0) = {-1,-2,-3,-4}
Conj(Q0) = {1,-2,-3,-4}
r + Q0 = {8,2,3,4}
Q0 + Q1 = {3,5,7,9}
Q1 + Q2 = {5,7,9,11}
Q1 * Q2 = {-56,16,24,26}
Q2 * Q1 = {-56,18,20,28}
Q1 * Q2 <> Q2 * Q1
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
irb(main):001:0> require 'quaternion' => true irb(main):002:0> q = Quaternion.new(1, 2, 3, 4) => (1, 2, 3, 4) irb(main):003:0> q1 = Quaternion.new(2, 3, 4, 5) => (2, 3, 4, 5) irb(main):004:0> q2 = Quaternion.new(3, 4, 5, 6) => (3, 4, 5, 6) irb(main):005:0> r = 7 => 7 irb(main):006:0> q.norm => 5.477225575051661 irb(main):007:0> q1.norm => 7.3484692283495345 irb(main):008:0> q2.norm => 9.273618495495704 irb(main):009:0> -q => (-1, -2, -3, -4) irb(main):010:0> q.conj => (1, -2, -3, -4) irb(main):011:0> q1 + q2 => (5, 7, 9, 11) irb(main):012:0> q2 + q1 => (5, 7, 9, 11) irb(main):013:0> q + r => (8, 2, 3, 4) irb(main):014:0> r + q => (8, 2, 3, 4) irb(main):015:0> q * r => (7, 14, 21, 28) irb(main):016:0> r * q => (7, 14, 21, 28) irb(main):017:0> q1 * q2 => (-56, 16, 24, 26) irb(main):018:0> q2 * q1 => (-56, 18, 20, 28) irb(main):019:0> q1 * q2 != q2 * q1 => true
Tcl
or
<lang tcl>package require TclOO
- 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]
}
}
- 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]"
- 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