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#
<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
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]
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 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
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
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