Quaternion type
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
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.
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>
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>
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
Haskell
<lang haskell>import Control.Monad import Control.Arrow import Data.List
data Quaternion = Q Double Double Double Double
deriving (Show, Ord, Eq)
realQ :: Quaternion -> Double realQ (Q r _ _ _) = r
imagQ :: Quaternion -> [Double] imagQ (Q _ i j k) = [i, j, k]
quaternionFromScalar s = Q s 0 0 0
listFromQ (Q a b c d) = [a,b,c,d] quaternionFromList [a, b, c, d] = Q a b c d
addQ, subQ, mulQ :: Quaternion -> Quaternion -> Quaternion addQ (Q a b c d) (Q p q r s) = Q (a+p) (b+q) (c+r) (d+s)
subQ (Q a b c d) (Q p q r s) = Q (a-p) (b-q) (c-r) (d-s)
mulQ (Q a b c d) (Q p q r s) =
Q (a*p - b*q - c*r - d*s) (a*q + b*p + c*s - d*r) (a*r - b*s + c*p + d*q) (a*s + b*r - c*q + d*p)
normQ = sqrt. sum. join (zipWith (*)). listFromQ
conjQ, negQ :: Quaternion -> Quaternion conjQ (Q a b c d) = Q a (-b) (-c) (-d)
negQ (Q a b c d) = Q (-a) (-b) (-c) (-d)</lang> To use with the Examples: <lang haskell>[q,q1,q2] = map quaternionFromList [[1..4],[2..5],[3..6]] -- a*b == b*a test :: Quaternion -> Quaternion -> Bool test a b = a `mulQ` b == b `mulQ` a</lang> Examples:
*Main> mulQ (Q 0 1 0 0) $ mulQ (Q 0 0 1 0) (Q 0 0 0 1) -- i*j*k Q (-1.0) 0.0 0.0 0.0 *Main> test q1 q2 False *Main> mulQ q1 q2 Q (-56.0) 16.0 24.0 26.0 *Main> flip mulQ q1 q2 Q (-56.0) 18.0 20.0 28.0 *Main> imagQ q [2.0,3.0,4.0]
J
Derived from the j wiki:
<lang j> NB. utilities
ip=: +/ .* NB. inner product t4=. (_1^#:0 10 9 12)*0 7 16 23 A.=i.4 toQ=: 4&{."1 :[: NB. real scalars -> quaternion
NB. task norm=: %:@ip~@toQ NB. | y neg=: -&toQ NB. - y and x - y conj=: 1 _1 _1 _1 * toQ NB. + y add=: +&toQ NB. x + y mul=: (ip t4 ip ])&toQ NB. x * y</lang>
Example use:
<lang> q=: 1 2 3 4
q1=: 2 3 4 5 q2=: 3 4 5 6 r=: 7 norm q
5.47723
neg q
_1 _2 _3 _4
conj q
1 _2 _3 _4
r add q
8 2 3 4
q1 add q2
5 7 9 11
r mul q
7 14 21 28
q1 mul q2
_56 16 24 26
q2 mul q1
_56 18 20 28</lang>
OCaml
The implementation as a file q.ml: <lang ocaml>type quaternion = float * float * float * float
let q a b c d = (a, b, c, d)
let to_real (r, _, _, _) = r let imag (_, i, j, k) = (i, j, k)
let quaternion_of_scalar s = (s, 0.0, 0.0, 0.0)
let to_list (a, b, c, d) = [a; b; c; d] let of_list = function [a; b; c; d] -> (a, b, c, d)
| _ -> invalid_arg "of_list"
let ( + ) = ( +. ) let ( - ) = ( -. ) let ( * ) = ( *. ) let ( / ) = ( /. )
let addr (a, b, c, d) r = (a+r, b, c, d) let mulr (a, b, c, d) r = (a*r, b*r, c*r, d*r)
let add (a, b, c, d) (p, q, r, s) = (a+p, b+q, c+r, d+s)
let sub (a, b, c, d) (p, q, r, s) = (a-p, b-q, c-r, d-s)
let mul (a, b, c, d) (p, q, r, s) =
( a*p - b*q - c*r - d*s, a*q + b*p + c*s - d*r, a*r - b*s + c*p + d*q, a*s + b*r - c*q + d*p )
let norm2 (a, b, c, d) =
( a * a + b * b + c * c + d * d )
let norm q = sqrt(norm2 q)
let conj (a, b, c, d) = (a, -. b, -. c, -. d) let neg (a, b, c, d) = (-. a, -. b, -. c, -. d)
let unit ((a, b, c, d) as q) =
let n = norm q in (a/n, b/n, c/n, d/n)
let reciprocal ((a, b, c, d) as q) =
let n2 = norm2 q in (a/n2, b/n2, c/n2, d/n2)</lang>
and the interface as a file q.mli: <lang ocaml>type quaternion = float * float * float * float
val q : float -> float -> float -> float -> quaternion val to_real : quaternion -> float val imag : quaternion -> float * float * float val quaternion_of_scalar : float -> quaternion val to_list : quaternion -> float list val of_list : float list -> quaternion val addr : quaternion -> float -> quaternion val mulr : quaternion -> float -> quaternion val add : quaternion -> quaternion -> quaternion val sub : quaternion -> quaternion -> quaternion val mul : quaternion -> quaternion -> quaternion val norm : quaternion -> float val conj : quaternion -> quaternion val neg : quaternion -> quaternion val unit : quaternion -> quaternion val reciprocal : quaternion -> quaternion</lang>
using this module in the interactive interpreter:
$ ocamlc -c q.mli $ ocamlc -c q.ml $ ocaml q.cmo Objective Caml version 3.11.2 # open Q ;; # let q1 = q 2.0 3.0 4.0 5.0 and q2 = q 3.0 4.0 5.0 6.0 ;; val q1 : Q.quaternion = (2., 3., 4., 5.) val q2 : Q.quaternion = (3., 4., 5., 6.) # (mul q1 q2) <> (mul q2 q1) ;; - : bool = true
Perl 6
<lang perl6>class Quaternion {
has Real ( $.r, $.i, $.j, $.k );
multi method new ( Real $r, Real $i, Real $j, Real $k ) { self.bless: *, :$r, :$i, :$j, :$k; } method Str ( ) { "$.r + {$.i}i + {$.j}j + {$.k}k" } method reals ( ) { $.r, $.i, $.j, $.k } method conj ( ) { self.new: $.r, -$.i, -$.j, -$.k } method norm ( ) { sqrt [+] self.reals X** 2 }
}
subset Qu of Quaternion; # Makes a short alias
multi sub infix:<eqv> ( Qu $a, Qu $b ) { [and] $a.reals Z== $b.reals } multi sub infix:<+> ( Qu $a, Real $b ) { $a.new: $b+$a.r, $a.i, $a.j, $a.k } multi sub infix:<+> ( Real $a, Qu $b ) { $b.new: $a+$b.r, $b.i, $b.j, $b.k } multi sub infix:<+> ( Qu $a, Qu $b ) { $a.new: |( $a.reals Z+ $b.reals ) } multi sub prefix:<-> ( Qu $a ) { $a.new: |( $a.reals X* -1 ) } multi sub infix:<*> ( Qu $a, Real $b ) { $a.new: |( $a.reals X* $b ) } multi sub infix:<*> ( Real $a, Qu $b ) { $b.new: |( $b.reals X* $a ) } multi sub infix:<*> ( Qu $a, Qu $b ) {
my @a_rijk = $a.reals; my ( $r, $i, $j, $k ) = $b.reals; return $a.new: ( [+] @a_rijk Z* $r, -$i, -$j, -$k ), # real ( [+] @a_rijk Z* $i, $r, $k, -$j ), # i ( [+] @a_rijk Z* $j, -$k, $r, $i ), # j ( [+] @a_rijk Z* $k, $j, -$i, $r ); # k
}
my Quaternion $q .= new: 1, 2, 3, 4; my Quaternion $q1 .= new: 2, 3, 4, 5; my Quaternion $q2 .= new: 3, 4, 5, 6; my $r = 7;
say "1) q norm = {$q.norm}"; say "2) -q = {-$q}"; say "3) q conj = {$q.conj}"; say "4) q + r = {$q + $r}"; say "5) q1 + q2 = {$q1 + $q2}"; say "6) q * r = {$q * $r}"; say "7) q1 * q2 = {$q1 * $q2}"; say "8) q1q2 { $q1 * $q2 eqv $q2 * $q1 ?? '==' !! '!=' } q2q1"; </lang>
Output:
1) q norm = 5.47722557505166 2) -q = -1 + -2i + -3j + -4k 3) q conj = 1 + -2i + -3j + -4k 4) q + r = 8 + 2i + 3j + 4k 5) q1 + q2 = 5 + 7i + 9j + 11k 6) q * r = 7 + 14i + 21j + 28k 7) q1 * q2 = -56 + 16i + 24j + 26k 8) q1q2 != q2q1
PureBasic
<lang PureBasic>Structure Quaternion
a.f: b.f c.f: d.f
EndStructure
Structure Quaternion2
Qa.Quaternion: Qb.Quaternion Qc.Quaternion: Qd.Quaternion
EndStructure
Procedure.f QNorm(*Q.Quaternion)
Protected Result.f Result=Sqr(Pow(*Q\a,2)+ Pow(*Q\b,2)+ Pow(*Q\c,2)+ Pow(*Q\d,2)) ProcedureReturn Result
EndProcedure
Procedure.i QNeg(*Q.Quaternion)
Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion)) If *B With *Q *B\a=-\a: *B\b=-\b *B\c=-\c: *B\d=-\d EndWith EndIf ProcedureReturn *B
EndProcedure
Procedure.i QConj(*Q.Quaternion)
Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion)) If *B With *Q *B\a= \a: *B\b=-\b *B\c=-\c: *B\d=-\d EndWith EndIf ProcedureReturn *B
EndProcedure
Procedure.i QAddReal(R.f, *Q.Quaternion)
Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion)) If *B With *Q *B\a=R+\a: *B\b= -\b *B\c= -\c: *B\d= -\d EndWith EndIf ProcedureReturn *B
EndProcedure
Procedure.i QAddQuaternion(*Q1.Quaternion, *Q2.Quaternion)
Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion)) If *B *B\a=*Q1\a + *Q2\a: *B\b=*Q1\b + *Q2\b *B\c=*Q1\c + *Q2\c: *B\d=*Q1\d + *Q2\d EndIf ProcedureReturn *B
EndProcedure
Procedure.i QMulReal_and_Quaternion(R.f, *Q.Quaternion)
Protected *B.Quaternion=AllocateMemory(SizeOf(Quaternion)) If *B *B\a=*Q\a * R: *B\b=*Q\b * R *B\c=*Q\c * R: *B\d=*Q\d * R EndIf ProcedureReturn *B
EndProcedure
Procedure.i QMulQuaternion(*Q1.Quaternion, *Q2.Quaternion)
Protected *B.Quaternion2=AllocateMemory(SizeOf(Quaternion2)) If *B With *B \Qa\a= *Q1\a * *Q2\a: \Qa\b=-*Q1\b * *Q2\b \Qa\c=-*Q1\c * *Q2\c: \Qa\d=-*Q1\d * *Q2\d ; \Qb\a= *Q1\a * *Q2\b: \Qb\b= *Q1\b * *Q2\a \Qb\c= *Q1\c * *Q2\d: \Qb\d=-*Q1\d * *Q2\c ; \Qc\a= *Q1\a * *Q2\c: \Qc\b=-*Q1\b * *Q2\d \Qc\c= *Q1\c * *Q2\a: \Qc\d= *Q1\d * *Q2\b ; \Qd\a= *Q1\a * *Q2\d: \Qd\b= *Q1\b * *Q2\c \Qd\c=-*Q1\c * *Q2\b: \Qd\d= *Q1\d * *Q2\a EndWith EndIf ProcedureReturn *B
EndProcedure</lang> Implementation & test <lang PureBasic>Macro Show(QQ,NN=0)
StrF(QQ\a,NN)+","+StrF(QQ\b,NN)+","+StrF(QQ\c,NN)+","+StrF(QQ\d,NN)
EndMacro
Define.Quaternion A, B, C Define.f r=7 Define *X.Quaternion, *Y.Quaternion2
A\a=1: A\b=2: A\c=3: A\d=4 B\a=2: B\b=3: B\c=4: B\d=5 C\a=3: C\b=4: C\c=5: C\d=6
Debug "Q1= {"+Show(A,0)+"}" Debug "Q2= {"+Show(B,0)+"}" Debug "Q3= {"+Show(C,0)+"}"
Debug "Normal of Q1= "+StrF(QNorm(@A))
- X=QNeg(@A)
Debug "Neg(Q1) ={"+Show(*X)+"}" : FreeMemory(*X)
- X=QConj(@A)
Debug "Conj(Q1) ={"+Show(*X)+"}" : FreeMemory(*X)
- X=QAddReal(r,@A)
Debug "r+A ={"+Show(*X)+"}" : FreeMemory(*X)
- X=QAddQuaternion(@A,@B)
Debug "Q1+Q2 ={"+Show(*X)+"}" : FreeMemory(*X)
- X=QAddQuaternion(@B,@C)
Debug "Q2+Q3 ={"+Show(*X)+"}" : FreeMemory(*X)
- Y=QMulQuaternion(@A,@B)
Debug "Q1*Q2 =" Debug "{{"+Show(*Y\Qa)+"}" Debug " {"+Show(*Y\Qb)+"}" Debug " {"+Show(*Y\Qc)+"}" Debug " {"+Show(*Y\Qd)+"}}" FreeMemory(*Y)
- Y=QMulQuaternion(@B,@A)
Debug "Q2*Q1 =" Debug "{{"+Show(*Y\Qa)+"}" Debug " {"+Show(*Y\Qb)+"}" Debug " {"+Show(*Y\Qc)+"}" Debug " {"+Show(*Y\Qd)+"}}" FreeMemory(*Y)</lang> Result
Q1= {1,2,3,4} Q2= {2,3,4,5} Q3= {3,4,5,6} Normal of Q1= 5.4772257805 Neg(Q1) ={-1,-2,-3,-4} Conj(Q1) ={1,-2,-3,-4} r+A ={8,-2,-3,-4} Q1+Q2 ={3,5,7,9} Q2+Q3 ={5,7,9,11} Q1*Q2 = {{2,-6,-12,-20} {3,4,15,-16} {4,-10,6,12} {5,8,-9,8}} Q2*Q1 = {{2,-6,-12,-20} {4,3,16,-15} {6,-12,4,10} {8,9,-8,5}}
Python
This example extends Pythons namedtuples to add extra functionality. <lang python>from collections import namedtuple import math
class Q(namedtuple('Quaternion', 'real, i, j, k')):
'Quaternion type: Q(real=0.0, i=0.0, j=0.0, k=0.0)'
__slots__ = ()
def __new__(_cls, real=0.0, i=0.0, j=0.0, k=0.0): 'Defaults all parts of quaternion to zero' return super().__new__(_cls, float(real), float(i), float(j), float(k))
def conjugate(self): return Q(self.real, -self.i, -self.j, -self.k)
def _norm2(self): return sum( x*x for x in self)
def norm(self): return math.sqrt(self._norm2())
def reciprocal(self): n2 = self._norm2() return Q(*(x / n2 for x in self.conjugate()))
def __str__(self): 'Shorter form of Quaternion as string' return 'Q(%g, %g, %g, %g)' % self
def __neg__(self): return Q(-self.real, -self.i, -self.j, -self.k)
def __add__(self, other): if type(other) == Q: return Q( *(s+o for s,o in zip(self, other)) ) try: f = float(other) except: return NotImplemented return Q(self.real + f, self.i, self.j, self.k)
def __radd__(self, other): return Q.__add__(self, other)
def __mul__(self, other): if type(other) == Q: a1,b1,c1,d1 = self a2,b2,c2,d2 = other return Q( a1*a2 - b1*b2 - c1*c2 - d1*d2, a1*b2 + b1*a2 + c1*d2 - d1*c2, a1*c2 - b1*d2 + c1*a2 + d1*b2, a1*d2 + b1*c2 - c1*b2 + d1*a2 ) try: f = float(other) except: return NotImplemented return Q(self.real * f, self.i * f, self.j * f, self.k * f)
def __rmul__(self, other): return Q.__mul__(self, other)
def __truediv__(self, other): if type(other) == Q: return self.__mul__(other.reciprocal()) try: f = float(other) except: return NotImplemented return Q(self.real / f, self.i / f, self.j / f, self.k / f)
def __rtruediv__(self, other): return other * self.reciprocal()
__div__, __rdiv__ = __truediv__, __rtruediv__
Quaternion = Q
q = Q(1, 2, 3, 4) q1 = Q(2, 3, 4, 5) q2 = Q(3, 4, 5, 6) r = 7</lang>
Continued shell session Run the above with the -i flag to python on the command line, or run with idle then continue in the shell as follows: <lang python>>>> q Quaternion(real=1.0, i=2.0, j=3.0, k=4.0) >>> q1 Quaternion(real=2.0, i=3.0, j=4.0, k=5.0) >>> q2 Quaternion(real=3.0, i=4.0, j=5.0, k=6.0) >>> r 7 >>> q.norm() 5.477225575051661 >>> q1.norm() 7.3484692283495345 >>> q2.norm() 9.273618495495704 >>> -q Quaternion(real=-1.0, i=-2.0, j=-3.0, k=-4.0) >>> q.conjugate() Quaternion(real=1.0, i=-2.0, j=-3.0, k=-4.0) >>> r + q Quaternion(real=8.0, i=2.0, j=3.0, k=4.0) >>> q + r Quaternion(real=8.0, i=2.0, j=3.0, k=4.0) >>> q1 + q2 Quaternion(real=5.0, i=7.0, j=9.0, k=11.0) >>> q2 + q1 Quaternion(real=5.0, i=7.0, j=9.0, k=11.0) >>> q * r Quaternion(real=7.0, i=14.0, j=21.0, k=28.0) >>> r * q Quaternion(real=7.0, i=14.0, j=21.0, k=28.0) >>> q1 * q2 Quaternion(real=-56.0, i=16.0, j=24.0, k=26.0) >>> q2 * q1 Quaternion(real=-56.0, i=18.0, j=20.0, k=28.0) >>> assert q1 * q2 != q2 * q1 >>> >>> i, j, k = Q(0,1,0,0), Q(0,0,1,0), Q(0,0,0,1) >>> i*i Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> j*j Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> k*k Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> i*j*k Quaternion(real=-1.0, i=0.0, j=0.0, k=0.0) >>> q1 / q2 Quaternion(real=0.7906976744186047, i=0.023255813953488358, j=-2.7755575615628914e-17, k=0.046511627906976744) >>> q1 / q2 * q2 Quaternion(real=2.0000000000000004, i=3.0000000000000004, j=4.000000000000001, k=5.000000000000001) >>> q2 * q1 / q2 Quaternion(real=2.0, i=3.465116279069768, j=3.906976744186047, k=4.767441860465116) >>> q1.reciprocal() * q1 Quaternion(real=0.9999999999999999, i=0.0, j=0.0, k=0.0) >>> q1 * q1.reciprocal() Quaternion(real=0.9999999999999999, i=0.0, j=0.0, k=0.0) >>> </lang>
Ruby
<lang ruby>class Quaternion
attr_reader :a, :b, :c, :d
def initialize a, b, c, d @a, @b, @c, @d = a, b, c, d end
def norm Math.sqrt(@a * @a + @b * @b + @c * @c + @d * @d) end
def -@ Quaternion.new(-@a, -@b, -@c, -@d) end
def conj Quaternion.new(@a, -@b, -@c, -@d) end
def + q if q.is_a? Quaternion Quaternion.new(@a + q.a, @b + q.b, @c + q.c, @d + q.d) elsif q.is_a? Numeric Quaternion.new(@a + q, @b, @c, @d) else q + self end end
def - q self + -q end
def * q if q.is_a? Quaternion Quaternion.new( @a * q.a - @b * q.b - @c * q.c - @d * q.d, @a * q.b - @b * q.a - @c * q.d - @d * q.c, @a * q.c - @b * q.d - @c * q.a - @d * q.b, @a * q.d - @b * q.c - @c * q.b - @d * q.a) elsif q.is_a? Numeric Quaternion.new(@a * q, @b * q, @c * q, @d * q) else q * self end end
def == q @a == q.a and @b == q.b and @c == q.c and @d == q.d end
def to_s "(#{@a}, #{@b}, #{@c}, #{@d})" end
# this is Ruby's way of handling commutativity def coerce other return self, other end
end</lang>
IRB session
>> require 'quaternions' => true >> q = Quaternion.new(1, 2, 3, 4) => #<Quaternion:0x7fcd33128a88 @b=2, @a=1, @d=4, @c=3> >> q1 = Quaternion.new(2, 3, 4, 5) => #<Quaternion:0x7fcd33122e80 @b=3, @a=2, @d=5, @c=4> >> q2 = Quaternion.new(3, 4, 5, 6) => #<Quaternion:0x7fcd3311d278 @b=4, @a=3, @d=6, @c=5> >> r = 7 => 7 >> q.norm => 5.47722557505166 >> q1.norm => 7.34846922834953 >> q2.norm => 9.2736184954957 >> -q => #<Quaternion:0x7fcd331149c0 @b=-2, @a=-1, @d=-4, @c=-3> >> q.conj => #<Quaternion:0x7fcd33112e18 @b=-2, @a=1, @d=-4, @c=-3> >> q1 + q2 => #<Quaternion:0x7fcd33110af0 @b=7, @a=5, @d=11, @c=9> >> q2 + q1 => #<Quaternion:0x7fcd3310e7c8 @b=7, @a=5, @d=11, @c=9> >> q + r => #<Quaternion:0x7fcd3310c590 @b=2, @a=8, @d=4, @c=3> >> r + q => #<Quaternion:0x7fcd3310a330 @b=2, @a=8, @d=4, @c=3> >> q * r => #<Quaternion:0x7fcd33108170 @b=14, @a=7, @d=28, @c=21> >> r * q => #<Quaternion:0x7fcd33105f88 @b=14, @a=7, @d=28, @c=21> >> q1 * q2 => #<Quaternion:0x7fcd33103cd8 @b=-50, @a=-56, @d=-34, @c=-40> >> q2 * q1 => #<Quaternion:0x7fcd33101a28 @b=-48, @a=-56, @d=-28, @c=-36> >> q1 * q2 != q2 * q1 => true
Tcl
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