# http://mathling.com/core/quaternion/vector  library module

http://mathling.com/core/quaternion/vector

Quaternions: vector implementation

Information on operations:
https://en.wikipedia.org/wiki/Quaternion
https://web.archive.org/web/20170705123142/http://www.lce.hut.fi/~ssarkka/pub/quat.pdf
https://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/functions/power/index.htm (and nearby)
https://math.stackexchange.com/questions/939229/unit-quaternion-to-a-scalar-power
https://math.stackexchange.com/questions/1499095/how-to-calculate-sin-cos-tan-of-a-quaternion
https://isidore.co/CalibreLibrary/Morais,%20Joao%20Pedro/Real%20Quaternionic%20Calculus%20Handbook%20(6130)/Real%20Quaternionic%20Calculus%20Handbook%20-%20Morais,%20Joao%20Pedro.pdf
https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles

November 2022
Status: Stable

### Imports

http://mathling.com/core/utilities
```import module namespace util="http://mathling.com/core/utilities"
at "../core/utilities.xqy"```
http://mathling.com/core/vector
```import module namespace v="http://mathling.com/core/vector"
at "../core/vector.xqy"```
http://mathling.com/geometric/affine
```import module namespace affine="http://mathling.com/geometric/affine"
at "../geo/affine.xqy"```
http://mathling.com/core/config
```import module namespace config="http://mathling.com/core/config"
at "../core/config.xqy"```
http://mathling.com/core/errors
```import module namespace errors="http://mathling.com/core/errors"
at "../core/errors.xqy"```

### Functions

#### ```Function: quaterniondeclare function quaternion(\$scalar as xs:double, \$i as xs:double, \$j as xs:double, \$k as xs:double) as xs:double*```

quaternion()
Construct a new quaternion

##### Params
• scalar as xs:double: the scalar part of the quaternion
• i as xs:double: i coordinate
• j as xs:double: j coordinate \$param \$k: k coordinate
• k as xs:double
##### Returns
• xs:double*: quaternion
```declare function this:quaternion(
\$scalar as xs:double,
\$i as xs:double,
\$j as xs:double,
\$k as xs:double
) as xs:double*
{
(\$scalar, \$i, \$j, \$k)
}```

#### ```Function: quaterniondeclare function quaternion(\$scalar as xs:double, \$vector as xs:double*) as xs:double*```

quaternion()
Construct a new quaternion using scalar + vector formulation

##### Params
• scalar as xs:double: the scalar part of the quaternion
• vector as xs:double*: the vector part of the quaternion
##### Returns
• xs:double*: quaternion
```declare function this:quaternion(
\$scalar as xs:double,
\$vector as xs:double*
) as xs:double*
{
(\$scalar, (\$vector[1],0)[1], (\$vector[2],0)[1], (\$vector[3],0)[1])
}```

#### `Function: quaterniondeclare function quaternion(\$coords as xs:double*) as xs:double*`

quaternion()
Construct a new quaternion using vector coordinates formulation

##### Params
• coords as xs:double*: coordinates of quaternion
##### Returns
• xs:double*: quaternion
```declare function this:quaternion(
\$coords as xs:double*
) as xs:double*
{
(\$coords)
}```

#### ```Function: quaterniondeclare function quaternion(\$roll-degrees as xs:double, \$pitch-degrees as xs:double, \$yaw-degrees as xs:double) as xs:double*```

quaternion()
Quaternion corresponding to given rotation angles

##### Params
• roll-degrees as xs:double degrees of roll
• pitch-degrees as xs:double degrees of pitch
• yaw-degrees as xs:double degrees of yaw
##### Returns
• xs:double*: quaternion
```declare function this:quaternion(
\$roll-degrees as xs:double,
\$pitch-degrees as xs:double,
\$yaw-degrees as xs:double
) as xs:double*
{
let \$cr := math:cos(util:radians(\$roll-degrees) div 2)
let \$sr := math:sin(util:radians(\$roll-degrees) div 2)
let \$cp := math:cos(util:radians(\$pitch-degrees) div 2)
let \$sp := math:sin(util:radians(\$pitch-degrees) div 2)
let \$cy := math:cos(util:radians(\$yaw-degrees) div 2)
let \$sy := math:sin(util:radians(\$yaw-degrees) div 2)
return (
this:quaternion(
\$cr * \$cp * \$cy + \$sr * \$sp * \$sy,
\$sr * \$cp * \$cy - \$cr * \$sp * \$sy,
\$cr * \$sp * \$cy + \$sr * \$cp * \$sy,
\$cr * \$cp * \$sy - \$sr * \$sp * \$cy
)
)
}```

#### ```Function: rotation-betweendeclare function rotation-between(\$from as xs:double*, \$to as xs:double*) as xs:double*```

rotation-between()
Compute the quaternion representing the rotation angle between two 3D
vectors.

##### Params
• from as xs:double*: starting vector
• to as xs:double*: ending vector
##### Returns
• xs:double*: quaternion representing the rotation
```declare function this:rotation-between(
\$from as xs:double*,
\$to as xs:double*
) as xs:double*
{
let \$dot := v:dot(\$from, \$to)
return
if (\$dot < -1.0 + \$config:ε) then (
let \$tv := v:cross( (1,0,0), \$from )
let \$tv :=
if (v:magnitude(\$tv) < \$config:ε)
then v:cross( (0,1,0), \$from )=>v:normalize()
else \$tv=>v:normalize()
return (
this:quaternion(0, \$tv) (: 180° :)
)
) else if (\$dot > 1.0 - \$config:ε) then (
\$this:one (: 0° :)
) else (
this:quaternion(1 + \$dot, v:cross(\$from, \$to))=>this:normalize()
)
}```

#### `Function: rotation3declare function rotation3(\$q as xs:double*) as xs:double*`

rotation3()
Construct the rotation matrix corresponding to the given quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: a matrix suitable for use by geom:apply-matrix3()
```declare function this:rotation3(
\$q as xs:double*
) as xs:double*
{
(: http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm :)
let \$qn := this:normalize(\$q)
let \$s := this:scalar(\$qn)
let \$v := this:vector(\$qn)
let \$qx := v:px(\$v)
let \$qy := v:py(\$v)
let \$qz := v:pz(\$v)
return (
1 - 2*\$qy*\$qy - 2*\$qz*\$qz, 2*\$qx*\$qy - 2*\$qz*\$s, 2*\$qx*\$qz + 2*\$qy*\$s, 0,
2*\$qx*\$qy + 2*\$qz*\$s, 1 - 2*\$qx*\$qx - 2*\$qz*\$qz, 2*\$qy*\$qz - 2*\$qx*\$s, 0,
2*\$qx*\$qz - 2*\$qy*\$s, 2*\$qy*\$qz + 2*\$qx*\$s, 1 - 2*\$qx*\$qx - 2*\$qy*\$qy, 0
)
}```

#### ```Function: rotation3declare function rotation3(\$qv as xs:double*, \$center as xs:double*) as xs:double*```

rotation3()
Construct the rotation matrix corresponding to the given quaternion

##### Params
• qv as xs:double*
• center as xs:double*: center of rotation, represented as a 3-vector
##### Returns
• xs:double*: a matrix suitable for use by geom:apply-matrix3()
```declare function this:rotation3(
\$qv as xs:double*,
\$center as xs:double*
) as xs:double*
{
affine:compose3(
affine:translation3(v:px(\$center), v:py(\$center), v:pz(\$center)),
affine:compose3(
this:rotation3(\$qv),
affine:translation3(-v:px(\$center), -v:py(\$center), -v:pz(\$center))
)
)
}```

#### `Function: to-anglesdeclare function to-angles(\$q as xs:double*) as xs:double*`

to-angles()
Convert quaternion to Euler angles.

##### Params
• q as xs:double*: quaternion
##### Returns
• xs:double*: roll, pitch, and yaw as degrees (same order as affine:rotation)
```declare function this:to-angles(
\$q as xs:double*
) as xs:double*
{
let \$qn := this:normalize(\$q)
let \$s := this:scalar(\$qn)
let \$v := this:vector(\$qn)
return (
(: roll x :)
let \$sinr_cosp := 2 * (\$s * v:px(\$v) + v:py(\$v) * v:pz(\$v))
let \$cosr_cosp := 1 - 2 * (v:px(\$v) * v:px(\$v) + v:py(\$v) * v:py(\$v))
return util:degrees(math:atan2(\$sinr_cosp, \$cosr_cosp))
,
(: pitch y :)
let \$sinp := math:sqrt(1 + 2 * (\$s * v:py(\$v) - v:px(\$v) * v:pz(\$v)))
let \$cosp := math:sqrt(1 - 2 * (\$s * v:py(\$v) - v:px(\$v) * v:pz(\$v)))
return util:degrees(2 * math:atan2(\$sinp, \$cosp) - math:pi() div 2)
,
(: yaw z :)
let \$siny_cosp := 2 * (\$s * v:pz(\$v) + v:px(\$v) * v:py(\$v))
let \$cosy_cosp := 1 - 2 * (v:py(\$v) * v:py(\$v) + v:pz(\$v) * v:pz(\$v))
return util:degrees(math:atan2(\$siny_cosp, \$cosy_cosp))
)
}```

#### `Function: as-quaterniondeclare function as-quaternion(\$scalar as xs:double) as xs:double*`

as-quaternion()
Cast a scalar number to a quaternion.

##### Params
• scalar as xs:double: the scalar value
##### Returns
• xs:double*: quaternion
```declare function this:as-quaternion(
\$scalar as xs:double
) as xs:double*
{
(\$scalar, 0, 0, 0)
}```

#### `Function: from-vectordeclare function from-vector(\$v as xs:double*) as xs:double*`

from-vector()
Convert a 3D vector into a rotation quaternion (normalized).

##### Params
• v as xs:double*: 3D vector
##### Returns
• xs:double*: normalized quaternion
```declare function this:from-vector(
\$v as xs:double*
) as xs:double*
{
let \$v := this:normalize(\$v)
let \$w := -math:sqrt(abs(1 - v:px(\$v)*v:px(\$v) - v:py(\$v)*v:py(\$v) - v:pz(\$v)*v:pz(\$v)))
return this:quaternion(\$w, \$v)
}```

#### `Function: is-scalardeclare function is-scalar(\$q as xs:double*) as xs:boolean`

is-scalar()
Test whether the quaternion has only a scalar part (Within ε).

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:boolean: true() if the quaternion only has a scalar part
```declare function this:is-scalar(
\$q as xs:double*
) as xs:boolean
{
every \$p in tail(\$q) satisfies abs(\$p) < \$config:ε
}```

#### `Function: as-scalardeclare function as-scalar(\$q as xs:double*) as xs:double`

as-scalar()
Cast a quaternion to a double, if possible. Raises error if not.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: double
```declare function this:as-scalar(
\$q as xs:double*
) as xs:double
{
if (this:is-scalar(\$q)) then this:ps(\$q)
else errors:error("UTIL-CAST", (\$q,'quaternion'))
}```
##### Errors

UTIL-CAST if the quaternion has non-zero i, j, or k parts

#### `Function: is-vectordeclare function is-vector(\$q as xs:double*) as xs:boolean`

is-vector()
Test whether the quaternion has only no scalar part (Within ε).

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:boolean: true() if there is no scalar part
```declare function this:is-vector(
\$q as xs:double*
) as xs:boolean
{
abs(this:ps(\$q)) < \$config:ε
}```

#### `Function: as-vectordeclare function as-vector(\$q as xs:double*) as xs:double*`

as-vector()
Cast a quaternion to its vector part, if possible. Raises error if not.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: the non-scalar part of the quaternion as a 3D vector (point)
```declare function this:as-vector(
\$q as xs:double*
) as xs:double*
{
if (this:is-vector(\$q)) then tail(\$q)
else errors:error("UTIL-CAST", (\$q,'quaternion'))
}```
##### Errors

UTIL-CAST if the quaternion has a scalar part

#### `Function: scalardeclare function scalar(\$q as xs:double*) as xs:double`

scalar()
Accessor for scalar part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: scalar part
```declare function this:scalar(\$q as xs:double*) as xs:double
{
v:px(\$q)
}```

#### `Function: vectordeclare function vector(\$q as xs:double*) as xs:double*`

vector()
Accessor for vector part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: vector part
```declare function this:vector(\$q as xs:double*) as xs:double*
{
v:py(\$q), v:pz(\$q), v:pw(\$q)
}```

#### `Function: angledeclare function angle(\$q as xs:double*) as xs:double`

angle()
cf arg() Returns as degrees
Quaternion (a,v) can be represented in polar form as ln(||q||) + e^(nθ)
cos(θ) = a/||q||, sin(θ) = ||v||/||q||
n = v/||v||
θ = acos(a/||q||) = asin(||v||/||q||)

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: phase angle in degrees
```declare function this:angle(\$q as xs:double*) as xs:double
{
if (this:scalar(\$q) = 0) then 90 else (
math:atan2(v:magnitude(this:vector(\$q)), this:scalar(\$q))=>
util:degrees()=>util:remap-degrees()
)
}```

#### `Function: argdeclare function arg(\$q as xs:double*) as xs:double`

arg()
Like angle(), but a name you see in math lit, and returns radians
aka phase() a name that speaks to me better

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: phase angle in radians
```declare function this:arg(\$q as xs:double*) as xs:double
{
this:phase(\$q)
}```

#### `Function: phasedeclare function phase(\$q as xs:double*) as xs:double`

phase()

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: phase angle in radians
```declare function this:phase(\$q as xs:double*) as xs:double
{
if (this:scalar(\$q) = 0) then math:pi() div 2
else math:atan(v:magnitude(this:vector(\$q)) div this:scalar(\$q))
}```

#### `Function: psdeclare function ps(\$q as xs:double*) as xs:double`

ps()
Accessor for real part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: raw real value
```declare function this:ps(\$q as xs:double*) as xs:double
{
v:px(\$q)
}```

#### `Function: pideclare function pi(\$q as xs:double*) as xs:double`

pi()
Accessor for i part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: raw i value
```declare function this:pi(\$q as xs:double*) as xs:double
{
v:py(\$q)
}```

#### `Function: pjdeclare function pj(\$q as xs:double*) as xs:double`

pj()
Accessor for j part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: raw j value
```declare function this:pj(\$q as xs:double*) as xs:double
{
v:pz(\$q)
}```

#### `Function: pkdeclare function pk(\$q as xs:double*) as xs:double`

pk()
Accessor for k part of the quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double: raw k value
```declare function this:pk(\$q as xs:double*) as xs:double
{
v:pw(\$q)
}```

#### `Function: sdeclare function s(\$q as xs:double*) as xs:integer`

s()
Accessor for real part of the quaternion, as integer.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:integer: raw real value, rounded
```declare function this:s(\$q as xs:double*) as xs:integer
{
v:x(\$q)
}```

#### `Function: ideclare function i(\$q as xs:double*) as xs:integer`

i()
Accessor for i part of the quaternion, as integer.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:integer: raw i value, rounded
```declare function this:i(\$q as xs:double*) as xs:integer
{
v:y(\$q)
}```

#### `Function: jdeclare function j(\$q as xs:double*) as xs:integer`

j()
Accessor for j part of the quaternion, as integer.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:integer: raw j value, rounded
```declare function this:j(\$q as xs:double*) as xs:integer
{
v:z(\$q)
}```

#### `Function: kdeclare function k(\$q as xs:double*) as xs:integer`

k()
Accessor for k part of the quaternion, as integer.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:integer: raw k value, rounded
```declare function this:k(\$q as xs:double*) as xs:integer
{
v:w(\$q)
}```

#### `Function: minusdeclare function minus(\$q as xs:double*) as xs:double*`

minus()
Negate the quaternion

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: -q
```declare function this:minus(
\$q as xs:double*
) as xs:double*
{
this:quaternion(\$q=>v:minus())
}```

#### ```Function: adddeclare function add(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double*```

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:double*: q1+q2
```declare function this:add(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
}```

#### ```Function: subdeclare function sub(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double*```

sub()
Subtract one quaternions from another.

##### Params
• q1 as xs:double*: the first quaternion
• q2 as xs:double*: the quaternion to subtract from the first quaternion
##### Returns
• xs:double*: q1-q2
```declare function this:sub(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
\$q1=>v:sub(\$q2)
}```

#### ```Function: multiplydeclare function multiply(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double*```

multiply()
Multiply two quaternions in the given order

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:double*: q1*q2
```declare function this:multiply(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
(: (a1 + b1i + c1j + d1k)*(a2 + b2i + c2j + d2k) =
(a1a2 - b1b2 - c1c2 - d1d2) +
(a1b2 + b1a2 + c1d2 - d1c2)i +
(a1c2 - b1d2 + c1a2 + d1b2)j +
(a1d2 + b1c2 - c1b2 + d1a2)k
:)
let \$a1 := this:ps(\$q1)
let \$b1 := this:pi(\$q1)
let \$c1 := this:pj(\$q1)
let \$d1 := this:pk(\$q1)
let \$a2 := this:ps(\$q2)
let \$b2 := this:pi(\$q2)
let \$c2 := this:pj(\$q2)
let \$d2 := this:pk(\$q2)
return (
this:quaternion(
(\$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)
)
)
}```

#### ```Function: dotdeclare function dot(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double```

dot()
Quaternion dot product

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:double: q1·q2
```declare function this:dot(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double
{
this:scalar(\$q1)*this:scalar(\$q2) + v:dot(this:vector(\$q1), this:vector(\$q2))
}```

#### ```Function: crossdeclare function cross(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double*```

cross()
Quaternion cross product

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:double*: q1Xq2
```declare function this:cross(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
this:quaternion(0, this:vector(\$q1=>this:multiply(\$q2)))
}```

#### ```Function: powdeclare function pow(\$q as xs:double*, \$pow as xs:double) as xs:double*```

pow()
Raise a quaternion to the given power.

##### Params
• q as xs:double*: the quaternion
• pow as xs:double: the power to raise it to
##### Returns
• xs:double*: q^pow
```declare function this:pow(
\$q as xs:double*,
\$pow as xs:double
) as xs:double*
{
if (this:is-scalar(\$q) and (this:ps(\$q) >= 0 or \$pow >= 1)) then (
this:as-quaternion(math:pow(this:ps(\$q), \$pow))
) else
switch (\$pow)
case 0 return \$this:one
case 1 return \$q
case 2 return
this:quaternion(
this:ps(\$q)*this:ps(\$q) - this:pi(\$q)*this:pi(\$q) - this:pj(\$q)*this:pj(\$q) - this:pk(\$q)*this:pk(\$q),
2*this:ps(\$q)*this:pi(\$q),
2*this:ps(\$q)*this:pj(\$q),
2*this:ps(\$q)*this:pk(\$q)
)
case 3 return
let \$a := this:ps(\$q)
let \$b := this:pi(\$q)
let \$c := this:pj(\$q)
let \$d := this:pk(\$q)
return (
this:quaternion(
math:pow(\$a,3) - 3*\$b*\$b*\$a - 3*\$c*\$c*\$a - 3*\$d*\$d*\$a,
3*\$a*\$a*\$b - math:pow(\$b,3) - 3*\$c*\$c*\$a - 3*\$d*\$d*\$a,
3*\$a*\$a*\$c - math:pow(\$c,3) - 3*\$b*\$b*\$c - 3*\$d*\$d*\$c,
3*\$a*\$a*\$d - math:pow(\$d,3) - 3*\$b*\$b*\$d - 3*\$c*\$c*\$d
)
)
default return (
this:exp(this:log(\$q)=>v:times(\$pow))
)
}```

#### `Function: logdeclare function log(\$q as xs:double*) as xs:double*`

log()
Natural log of quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: log(q)
```declare function this:log(
\$q as xs:double*
) as xs:double*
{
let \$m := this:magnitude(\$q)
return (
this:quaternion(
math:log(\$m),
(this:vector(\$q)=>v:normalize())=>
v:times(math:acos(this:scalar(\$q) div \$m))
)
)
}```

#### `Function: expdeclare function exp(\$q as xs:double*) as xs:double*`

exp()
Exponential of q

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: e^q
```declare function this:exp(
\$q as xs:double*
) as xs:double*
{
let \$ea := math:exp(this:scalar(\$q))
let \$nv := v:magnitude(this:vector(\$q))
return (
this:quaternion(
\$ea * math:cos(\$nv),
(this:vector(\$q)=>v:normalize())=>v:times(\$ea * math:sin(\$nv))
)
)
}```

#### `Function: sqrtdeclare function sqrt(\$q as xs:double*) as xs:double*`

sqrt()
Principal square root of quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: √q
```declare function this:sqrt(
\$q as xs:double*
) as xs:double*
{
if (this:is-scalar(\$q)) then (
if (this:ps(\$q) >= 0) then this:as-quaternion(math:sqrt(this:ps(\$q)))
else this:quaternion(0, math:sqrt(abs(this:ps(\$q))), 0, 0)
) else (
this:quaternion(
math:sqrt((this:magnitude(\$q) + this:ps(\$q)) div 2),
v:normalize(this:vector(\$q))=>
v:times(math:sqrt((this:magnitude(\$q) - this:ps(\$q)) div 2))
)
)
}```

#### `Function: sindeclare function sin(\$q as xs:double*) as xs:double*`

sin()
Sine of quaternion. Warning: some tests involving q:sin() are still failing
(although many are not). Still working it out.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: sin(q)
```declare function this:sin(
\$q as xs:double*
) as xs:double*
{
(:
: sin(q) = -1/2 n(e^(q n) - e^(-q n)) |v(q)|!=0; sin(s) otherwise
: where n = v(q)/|v(q)| i.e. normalized vector(q)
: which is -n sinh(qn)
:)
if (abs(v:magnitude(this:vector(\$q))) > \$config:ε) then (
let \$nv := this:quaternion(0, this:vector(\$q)=>v:normalize())
return this:minus(\$nv)=>this:multiply( this:sinh(\$q=>this:multiply(\$nv)) )
) else (
this:as-quaternion(math:sin(this:scalar(\$q)))
)
}```

#### `Function: cosdeclare function cos(\$q as xs:double*) as xs:double*`

cos()
Cosine of quaternion.

##### Params
• q as xs:double*: the quaternion
##### Returns
• xs:double*: cos(q)
```declare function this:cos(
\$q as xs:double*
) as xs:double*
{
(:
: cos(q) = 1/2 (e^(q n) + e^(-q n)) |q|!=0; cos(s) otherwise
: where n = v(q)/|v(q)| i.e. normalized vector(q)
: which is cosh(qn)
:)
if (abs(v:magnitude(this:vector(\$q))) > \$config:ε) then (
let \$nv := this:quaternion(0, this:vector(\$q)=>v:normalize())
return this:cosh( \$q=>this:multiply(\$nv) )
) else (
this:as-quaternion(math:cos(this:scalar(\$q)))
)
}```

#### `Function: sinhdeclare function sinh(\$x as xs:double*) as xs:double*`

sinh()
Hyperbolic sine.

##### Params
• x as xs:double*
##### Returns
• xs:double*: sinh(q)
```declare function this:sinh(\$x as xs:double*) as xs:double*
{
(this:exp(\$x)=>this:sub( this:exp(this:minus(\$x)) ))=>this:times(1 div 2.0)
}```

#### `Function: coshdeclare function cosh(\$x as xs:double*) as xs:double*`

cosh()
Hyperbolic cosine.

##### Params
• x as xs:double*
##### Returns
• xs:double*: cosh(q)
```declare function this:cosh(\$x as xs:double*) as xs:double*
{
}```

#### ```Function: timesdeclare function times(\$q as xs:double*, \$k as xs:double) as xs:double*```

times()
Multiply a quaternion by a constant.

##### Params
• q as xs:double*: the quaternion
• k as xs:double: the constant to multiply by
##### Returns
• xs:double*: k*q
```declare function this:times(
\$q as xs:double*,
\$k as xs:double
) as xs:double*
{
\$q=>v:times(\$k)
}```

#### ```Function: plusdeclare function plus(\$q as xs:double*, \$k as xs:double) as xs:double*```

plus()
Add a constant to a quaternion. (i.e. add quaternion (\$k, 0, 0, 0))

##### Params
• q as xs:double*: the quaternion
• k as xs:double: the constant to add
##### Returns
• xs:double*: q+k
```declare function this:plus(
\$q as xs:double*,
\$k as xs:double
) as xs:double*
{
}```

#### ```Function: dividedeclare function divide(\$q1 as xs:double*, \$q2 as xs:double*) as xs:double*```

divide()
Divide two quaternions q1/q2 as q1*q2^-1

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:double*: q1/q2
```declare function this:divide(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
\$q1=>this:multiply(this:reciprocal(\$q2))
}```

#### `Function: magnitudedeclare function magnitude(\$q as xs:double*) as xs:double`

magnitude()
Magnitude of a quaternion; aka modulus()

##### Params
• q as xs:double*: the quaternions
##### Returns
• xs:double: length of quaternion vector
```declare function this:magnitude(
\$q as xs:double*
) as xs:double
{
v:magnitude(\$q)
}```

#### `Function: modulusdeclare function modulus(\$q as xs:double*) as xs:double`

modulus()
Same as magnitude(), name you see in math papers

##### Params
• q as xs:double*: the quaternions
##### Returns
• xs:double: modulus
```declare function this:modulus(
\$q as xs:double*
) as xs:double
{
v:magnitude(\$q)
}```

#### `Function: normalizedeclare function normalize(\$q as xs:double*) as xs:double*`

normalize()
Normalize the quaternion so magnitude is 1. aka sgn()
Fails for (0,0,0,0)

##### Params
• q as xs:double*: the quaternion to normalize
##### Returns
• xs:double*: s normalized quaternion
```declare function this:normalize(\$q as xs:double*) as xs:double*
{
let \$l := this:magnitude(\$q)
return (
this:quaternion(
this:ps(\$q) div \$l,
this:pi(\$q) div \$l,
this:pj(\$q) div \$l,
this:pk(\$q) div \$l
)
)
}```
##### Errors

Divide by zero for the origin

#### `Function: sgndeclare function sgn(\$q as xs:double*) as xs:double*`

sgn()
Same as normalize(); name you see in math papers
Doesn't speak to me.

##### Params
• q as xs:double*: the quaternion to normalize
##### Returns
• xs:double*: s normalized quaternion
```declare function this:sgn(\$q as xs:double*) as xs:double*
{
this:normalize(\$q)
}```
##### Errors

Divide by zero for the origin

#### `Function: conjugatedeclare function conjugate(\$q as xs:double*) as xs:double*`

conjugate()
Return the quaternion conjugate

##### Params
• q as xs:double*: quaternion
##### Returns
• xs:double*: conj(q)
```declare function this:conjugate(\$q as xs:double*) as xs:double*
{
this:quaternion(this:ps(\$q), -this:pi(\$q), -this:pj(\$q), -this:pk(\$q))
}```

#### `Function: reciprocaldeclare function reciprocal(\$q as xs:double*) as xs:double*`

reciprocal()
The reciprocal of q: q^-1 = conj(q)/||q||²

##### Params
• q as xs:double*: quaternion
##### Returns
• xs:double*: q^-1
```declare function this:reciprocal(\$q as xs:double*) as xs:double*
{
(: Avoid the extra sqrt :)
let \$n2 :=
this:ps(\$q)*this:ps(\$q) +
this:pi(\$q)*this:pi(\$q) +
this:pj(\$q)*this:pj(\$q) +
this:pk(\$q)*this:pk(\$q)
return this:conjugate(\$q)=>this:times(1 div \$n2)
}```

#### `Function: absdeclare function abs(\$q as xs:double*) as xs:double*`

abs()
The absolute value of z

##### Params
• q as xs:double*: quaternion
##### Returns
• xs:double*: abs(q)
```declare function this:abs(\$q as xs:double*) as xs:double*
{
this:quaternion(\$q=>v:abs())
}```

#### `Function: quotedeclare function quote(\$qs as xs:double*) as xs:string`

quote()
Return a string value of the quaternion

##### Params
• qs as xs:double*
##### Returns
• xs:string: string
```declare function this:quote(\$qs as xs:double*) as xs:string
{
string-join(
let \$chunks :=
for \$q at \$i in \$qs return (
if (\$i mod 4 = 0) then (
array {\$qs[\$i - 3], \$qs[\$i - 2], \$qs[\$i - 1], \$q}
) else ()
)
for \$chunk in \$chunks
let \$q := \$chunk?*
return (
let \$s := this:scalar(\$q)
let \$v := this:vector(\$q)
let \$ps := if (\$s != 0 or (every \$c in \$v satisfies \$c = 0)) then string(\$s) else ()
let \$pi :=
if (\$v[1] = 0) then ()
else if (\$v[1] = -1) then "-i"
else if (empty(\$ps)) then (if (\$v[1]=1) then "i" else \$v[1]||"i")
else if (\$v[1] < 0) then \$v[1]||"i"
else if (\$v[1] = 1) then "+i"
else "+"||\$v[1]||"i"
let \$pj :=
if (\$v[2] = 0) then ()
else if (\$v[2] = -1) then "-j"
else if (empty((\$ps, \$pi))) then (if (\$v[2]=1) then "j" else \$v[2]||"j")
else if (\$v[2] < 0) then \$v[2]||"j"
else if (\$v[2] = 1) then "+j"
else "+"||\$v[2]||"j"
let \$pk :=
if (\$v[3] = 0) then ()
else if (\$v[3] = -1) then "-k"
else if (empty((\$ps, \$pi, \$pj))) then (if (\$v[3]=1) then "k" else \$v[3]||"k")
else if (\$v[3] < 0) then \$v[3]||"k"
else if (\$v[3] = 1) then "+k"
else "+"||\$v[3]||"k"
return (
if (\$s=0 and exists((\$pi,\$pj,\$pk)))
then \$pi||\$pj||\$pk
else \$ps||\$pi||\$pj||\$pk
)
)
,
" "
)
}```

#### `Function: unquotedeclare function unquote(\$qs as xs:string) as xs:double*`

unquote()
Parse a string value of a quaternion to produce that quaternion.
Number is of the form 33+12i+6j+8k

##### Params
• qs as xs:string: the string representation of the number
##### Returns
• xs:double*: quaternion
```declare function this:unquote(\$qs as xs:string) as xs:double*
{
let \$analysis :=
analyze-string(
replace(\$qs, " ", ""),
"([+-]?((\d+([.]\d)?[ijk]?)|(\d*[.]\d+[ijk?])|([ijk])))"
)
let \$parts := \$analysis/fn:match
return (
if (count(\$parts) > 4 or exists(\$analysis/fn:non-match))
else (
let \$six :=
(for \$part at \$i in \$parts where not(matches(\$part, "[ijk]")) return \$i)
let \$iix :=
for \$part at \$i in \$parts where ends-with(\$part, "i") return \$i
let \$jix :=
for \$part at \$i in \$parts where ends-with(\$part, "j") return \$i
let \$kix :=
for \$part at \$i in \$parts where ends-with(\$part, "k") return \$i
let \$number :=
function(\$ix as xs:integer?, \$v as xs:string) as xs:double {
if (empty(\$ix)) then 0 else (
let \$s := string(\$parts[\$ix])
let \$prefix := substring-before(\$s, \$v)
return (
if (empty(\$s)) then 0
else if (\$prefix=("","+")) then 1
else if (\$prefix="-") then -1
else number(\$prefix)
)
)
}
return (
if (count(\$six) > 1 or count(\$iix) > 1 or count(\$jix) > 1 or count(\$kix) > 1)
else (
this:quaternion(
if (exists(\$six)) then number(string(\$parts[\$six])) else 0,
\$number(\$iix, "i"),
\$number(\$jix, "j"),
\$number(\$kix, "k")
)
)
)
)
)
}```

#### `Function: samedeclare function same(\$q1 as xs:double*, \$q2 as xs:double*) as xs:boolean`

same()
Are the two quaternions the same?

##### Params
• q1 as xs:double*: one quaternion
• q2 as xs:double*: other quaternion
##### Returns
• xs:boolean: true() if the values are the same
```declare function this:same(\$q1 as xs:double*, \$q2 as xs:double*) as xs:boolean
{
this:ps(\$q1)=this:ps(\$q2) and
this:pi(\$q1)=this:pi(\$q2) and
this:pj(\$q1)=this:pj(\$q2) and
this:pk(\$q1)=this:pk(\$q2)
}```

#### `Function: decimaldeclare function decimal(\$qs as xs:double*, \$digits as xs:integer) as xs:double*`

decimal()
Cast the values to a values with the given number of digits.

##### Params
• qs as xs:double*: the quaternions
• digits as xs:integer: how many digits to preserve
##### Returns
• xs:double*: complex numbers
```declare function this:decimal(\$qs as xs:double*, \$digits as xs:integer) as xs:double*
{
v:decimal(\$qs, \$digits)
}```

### Original Source Code

```xquery version "3.1";
(:~
: Quaternions: vector implementation
:
: Information on operations:
: https://en.wikipedia.org/wiki/Quaternion
: https://web.archive.org/web/20170705123142/http://www.lce.hut.fi/~ssarkka/pub/quat.pdf
: https://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions/functions/power/index.htm (and nearby)
: https://math.stackexchange.com/questions/939229/unit-quaternion-to-a-scalar-power
: https://math.stackexchange.com/questions/1499095/how-to-calculate-sin-cos-tan-of-a-quaternion
: https://isidore.co/CalibreLibrary/Morais,%20Joao%20Pedro/Real%20Quaternionic%20Calculus%20Handbook%20(6130)/Real%20Quaternionic%20Calculus%20Handbook%20-%20Morais,%20Joao%20Pedro.pdf
: https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
:
: @since November 2022
: @custom:Status Stable
:)
module namespace this="http://mathling.com/core/quaternion/vector";

import module namespace config="http://mathling.com/core/config"
at "../core/config.xqy";
import module namespace errors="http://mathling.com/core/errors"
at "../core/errors.xqy";
import module namespace util="http://mathling.com/core/utilities"
at "../core/utilities.xqy";
import module namespace v="http://mathling.com/core/vector"
at "../core/vector.xqy";
import module namespace affine="http://mathling.com/geometric/affine"
at "../geo/affine.xqy";

declare namespace map="http://www.w3.org/2005/xpath-functions/map";
declare namespace math="http://www.w3.org/2005/xpath-functions/math";

declare variable \$this:i as xs:double* := this:quaternion(0, 1, 0, 0);
declare variable \$this:j as xs:double* := this:quaternion(0, 0, 1, 0);
declare variable \$this:k as xs:double* := this:quaternion(0, 0, 0, 1);
declare variable \$this:one as xs:double* := this:quaternion(1, 0, 0, 0);

(:~
: quaternion()
: Construct a new quaternion
:
: @param \$scalar: the scalar part of the quaternion
: @param \$i: i coordinate
: @param \$j: j coordinate
: \$param \$k: k coordinate
: @return quaternion
:)
declare function this:quaternion(
\$scalar as xs:double,
\$i as xs:double,
\$j as xs:double,
\$k as xs:double
) as xs:double*
{
(\$scalar, \$i, \$j, \$k)
};

(:~
: quaternion()
: Construct a new quaternion using scalar + vector formulation
:
: @param \$scalar: the scalar part of the quaternion
: @param \$vector: the vector part of the quaternion
: @return quaternion
:)
declare function this:quaternion(
\$scalar as xs:double,
\$vector as xs:double*
) as xs:double*
{
(\$scalar, (\$vector[1],0)[1], (\$vector[2],0)[1], (\$vector[3],0)[1])
};

(:~
: quaternion()
: Construct a new quaternion using vector coordinates formulation
:
: @param \$coords: coordinates of quaternion
: @return quaternion
:)
declare function this:quaternion(
\$coords as xs:double*
) as xs:double*
{
(\$coords)
};

(:~
: quaternion()
: Quaternion corresponding to given rotation angles
:
: @param \$roll-degrees degrees of roll
: @param \$pitch-degrees degrees of pitch
: @param \$yaw-degrees degrees of yaw
: @return quaternion
:)
declare function this:quaternion(
\$roll-degrees as xs:double,
\$pitch-degrees as xs:double,
\$yaw-degrees as xs:double
) as xs:double*
{
let \$cr := math:cos(util:radians(\$roll-degrees) div 2)
let \$sr := math:sin(util:radians(\$roll-degrees) div 2)
let \$cp := math:cos(util:radians(\$pitch-degrees) div 2)
let \$sp := math:sin(util:radians(\$pitch-degrees) div 2)
let \$cy := math:cos(util:radians(\$yaw-degrees) div 2)
let \$sy := math:sin(util:radians(\$yaw-degrees) div 2)
return (
this:quaternion(
\$cr * \$cp * \$cy + \$sr * \$sp * \$sy,
\$sr * \$cp * \$cy - \$cr * \$sp * \$sy,
\$cr * \$sp * \$cy + \$sr * \$cp * \$sy,
\$cr * \$cp * \$sy - \$sr * \$sp * \$cy
)
)
};

(:~
: rotation-between()
: Compute the quaternion representing the rotation angle between two 3D
: vectors.
: @param \$from: starting vector
: @param \$to: ending vector
: @return quaternion representing the rotation
:)
declare function this:rotation-between(
\$from as xs:double*,
\$to as xs:double*
) as xs:double*
{
let \$dot := v:dot(\$from, \$to)
return
if (\$dot < -1.0 + \$config:ε) then (
let \$tv := v:cross( (1,0,0), \$from )
let \$tv :=
if (v:magnitude(\$tv) < \$config:ε)
then v:cross( (0,1,0), \$from )=>v:normalize()
else \$tv=>v:normalize()
return (
this:quaternion(0, \$tv) (: 180° :)
)
) else if (\$dot > 1.0 - \$config:ε) then (
\$this:one (: 0° :)
) else (
this:quaternion(1 + \$dot, v:cross(\$from, \$to))=>this:normalize()
)
};

(:~
: rotation3()
: Construct the rotation matrix corresponding to the given quaternion.
: @param \$q: the quaternion
: @return a matrix suitable for use by geom:apply-matrix3()
:)
declare function this:rotation3(
\$q as xs:double*
) as xs:double*
{
(: http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToMatrix/index.htm :)
let \$qn := this:normalize(\$q)
let \$s := this:scalar(\$qn)
let \$v := this:vector(\$qn)
let \$qx := v:px(\$v)
let \$qy := v:py(\$v)
let \$qz := v:pz(\$v)
return (
1 - 2*\$qy*\$qy - 2*\$qz*\$qz, 2*\$qx*\$qy - 2*\$qz*\$s, 2*\$qx*\$qz + 2*\$qy*\$s, 0,
2*\$qx*\$qy + 2*\$qz*\$s, 1 - 2*\$qx*\$qx - 2*\$qz*\$qz, 2*\$qy*\$qz - 2*\$qx*\$s, 0,
2*\$qx*\$qz - 2*\$qy*\$s, 2*\$qy*\$qz + 2*\$qx*\$s, 1 - 2*\$qx*\$qx - 2*\$qy*\$qy, 0
)
};

(:~
: rotation3()
: Construct the rotation matrix corresponding to the given quaternion
: rotated about the given point.
: @param \$q: the quaternion
: @param \$center: center of rotation, represented as a 3-vector
: @return a matrix suitable for use by geom:apply-matrix3()
:)
declare function this:rotation3(
\$qv as xs:double*,
\$center as xs:double*
) as xs:double*
{
affine:compose3(
affine:translation3(v:px(\$center), v:py(\$center), v:pz(\$center)),
affine:compose3(
this:rotation3(\$qv),
affine:translation3(-v:px(\$center), -v:py(\$center), -v:pz(\$center))
)
)
};

(:~
: to-angles()
: Convert quaternion to Euler angles.
:
: @param \$q: quaternion
: @return roll, pitch, and yaw as degrees (same order as affine:rotation)
:)
declare function this:to-angles(
\$q as xs:double*
) as xs:double*
{
let \$qn := this:normalize(\$q)
let \$s := this:scalar(\$qn)
let \$v := this:vector(\$qn)
return (
(: roll x :)
let \$sinr_cosp := 2 * (\$s * v:px(\$v) + v:py(\$v) * v:pz(\$v))
let \$cosr_cosp := 1 - 2 * (v:px(\$v) * v:px(\$v) + v:py(\$v) * v:py(\$v))
return util:degrees(math:atan2(\$sinr_cosp, \$cosr_cosp))
,
(: pitch y :)
let \$sinp := math:sqrt(1 + 2 * (\$s * v:py(\$v) - v:px(\$v) * v:pz(\$v)))
let \$cosp := math:sqrt(1 - 2 * (\$s * v:py(\$v) - v:px(\$v) * v:pz(\$v)))
return util:degrees(2 * math:atan2(\$sinp, \$cosp) - math:pi() div 2)
,
(: yaw z :)
let \$siny_cosp := 2 * (\$s * v:pz(\$v) + v:px(\$v) * v:py(\$v))
let \$cosy_cosp := 1 - 2 * (v:py(\$v) * v:py(\$v) + v:pz(\$v) * v:pz(\$v))
return util:degrees(math:atan2(\$siny_cosp, \$cosy_cosp))
)
};

(:~
: as-quaternion()
: Cast a scalar number to a quaternion.
:
: @param \$scalar: the scalar value
: @return quaternion
:)
declare function this:as-quaternion(
\$scalar as xs:double
) as xs:double*
{
(\$scalar, 0, 0, 0)
};

(:~
: from-vector()
: Convert a 3D vector into a rotation quaternion (normalized).
: @param \$v: 3D vector
: @return normalized quaternion
:)
declare function this:from-vector(
\$v as xs:double*
) as xs:double*
{
let \$v := this:normalize(\$v)
let \$w := -math:sqrt(abs(1 - v:px(\$v)*v:px(\$v) - v:py(\$v)*v:py(\$v) - v:pz(\$v)*v:pz(\$v)))
return this:quaternion(\$w, \$v)
};

(:~
: is-scalar()
: Test whether the quaternion has only a scalar part (Within ε).
:
: @param \$q: the quaternion
: @return true() if the quaternion only has a scalar part
:)
declare function this:is-scalar(
\$q as xs:double*
) as xs:boolean
{
every \$p in tail(\$q) satisfies abs(\$p) < \$config:ε
};

(:~
: as-scalar()
: Cast a quaternion to a double, if possible. Raises error if not.
:
: @param \$q: the quaternion
: @return double
: @error UTIL-CAST if the quaternion has non-zero i, j, or k parts
:)
declare function this:as-scalar(
\$q as xs:double*
) as xs:double
{
if (this:is-scalar(\$q)) then this:ps(\$q)
else errors:error("UTIL-CAST", (\$q,'quaternion'))
};

(:~
: is-vector()
: Test whether the quaternion has only no scalar part (Within ε).
:
: @param \$q: the quaternion
: @return true() if there is no scalar part
:)
declare function this:is-vector(
\$q as xs:double*
) as xs:boolean
{
abs(this:ps(\$q)) < \$config:ε
};

(:~
: as-vector()
: Cast a quaternion to its vector part, if possible. Raises error if not.
:
: @param \$q: the quaternion
: @return the non-scalar part of the quaternion as a 3D vector (point)
: @error UTIL-CAST if the quaternion has a scalar part
:)
declare function this:as-vector(
\$q as xs:double*
) as xs:double*
{
if (this:is-vector(\$q)) then tail(\$q)
else errors:error("UTIL-CAST", (\$q,'quaternion'))
};

(:~
: scalar()
: Accessor for scalar part of the quaternion.
:
: @param \$q: the quaternion
: @return scalar part
:)
declare function this:scalar(\$q as xs:double*) as xs:double
{
v:px(\$q)
};

(:~
: vector()
: Accessor for vector part of the quaternion.
:
: @param \$q: the quaternion
: @return vector part
:)
declare function this:vector(\$q as xs:double*) as xs:double*
{
v:py(\$q), v:pz(\$q), v:pw(\$q)
};

(:~
: angle()
: cf arg() Returns as degrees
: Quaternion (a,v) can be represented in polar form as ln(||q||) + e^(nθ)
: cos(θ) = a/||q||, sin(θ) = ||v||/||q||
: n = v/||v||
: θ = acos(a/||q||) = asin(||v||/||q||)
:
: @param \$q: the quaternion
: @return phase angle in degrees
:)
declare function this:angle(\$q as xs:double*) as xs:double
{
if (this:scalar(\$q) = 0) then 90 else (
math:atan2(v:magnitude(this:vector(\$q)), this:scalar(\$q))=>
util:degrees()=>util:remap-degrees()
)
};

(:~
: arg()
: Like angle(), but a name you see in math lit, and returns radians
: aka phase() a name that speaks to me better
:
: @param \$q: the quaternion
: @return phase angle in radians
:)
declare function this:arg(\$q as xs:double*) as xs:double
{
this:phase(\$q)
};

(:~
: phase()
: Like angle(), but returns radians.
:
: @param \$q: the quaternion
: @return phase angle in radians
:)
declare function this:phase(\$q as xs:double*) as xs:double
{
if (this:scalar(\$q) = 0) then math:pi() div 2
else math:atan(v:magnitude(this:vector(\$q)) div this:scalar(\$q))
};

(:~
: ps()
: Accessor for real part of the quaternion.
:
: @param \$q: the quaternion
: @return raw real value
:)
declare function this:ps(\$q as xs:double*) as xs:double
{
v:px(\$q)
};

(:~
: pi()
: Accessor for i part of the quaternion.
:
: @param \$q: the quaternion
: @return raw i value
:)
declare function this:pi(\$q as xs:double*) as xs:double
{
v:py(\$q)
};

(:~
: pj()
: Accessor for j part of the quaternion.
:
: @param \$q: the quaternion
: @return raw j value
:)
declare function this:pj(\$q as xs:double*) as xs:double
{
v:pz(\$q)
};

(:~
: pk()
: Accessor for k part of the quaternion.
:
: @param \$q: the quaternion
: @return raw k value
:)
declare function this:pk(\$q as xs:double*) as xs:double
{
v:pw(\$q)
};

(:~
: s()
: Accessor for real part of the quaternion, as integer.
:
: @param \$q: the quaternion
: @return raw real value, rounded
:)
declare function this:s(\$q as xs:double*) as xs:integer
{
v:x(\$q)
};

(:~
: i()
: Accessor for i part of the quaternion, as integer.
:
: @param \$q: the quaternion
: @return raw i value, rounded
:)
declare function this:i(\$q as xs:double*) as xs:integer
{
v:y(\$q)
};

(:~
: j()
: Accessor for j part of the quaternion, as integer.
:
: @param \$q: the quaternion
: @return raw j value, rounded
:)
declare function this:j(\$q as xs:double*) as xs:integer
{
v:z(\$q)
};

(:~
: k()
: Accessor for k part of the quaternion, as integer.
:
: @param \$q: the quaternion
: @return raw k value, rounded
:)
declare function this:k(\$q as xs:double*) as xs:integer
{
v:w(\$q)
};

(:~
: minus()
: Negate the quaternion
:
: @param \$q: the quaternion
: @return -q
:)
declare function this:minus(
\$q as xs:double*
) as xs:double*
{
this:quaternion(\$q=>v:minus())
};

(:~
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return q1+q2
:)
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
};

(:~
: sub()
: Subtract one quaternions from another.
:
: @param \$q1: the first quaternion
: @param \$q2: the quaternion to subtract from the first quaternion
: @return q1-q2
:)
declare function this:sub(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
\$q1=>v:sub(\$q2)
};

(:~
: multiply()
: Multiply two quaternions in the given order
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return q1*q2
:)
declare function this:multiply(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
(: (a1 + b1i + c1j + d1k)*(a2 + b2i + c2j + d2k) =
(a1a2 - b1b2 - c1c2 - d1d2) +
(a1b2 + b1a2 + c1d2 - d1c2)i +
(a1c2 - b1d2 + c1a2 + d1b2)j +
(a1d2 + b1c2 - c1b2 + d1a2)k
:)
let \$a1 := this:ps(\$q1)
let \$b1 := this:pi(\$q1)
let \$c1 := this:pj(\$q1)
let \$d1 := this:pk(\$q1)
let \$a2 := this:ps(\$q2)
let \$b2 := this:pi(\$q2)
let \$c2 := this:pj(\$q2)
let \$d2 := this:pk(\$q2)
return (
this:quaternion(
(\$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)
)
)
};

(:~
: dot()
: Quaternion dot product
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return q1·q2
:)
declare function this:dot(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double
{
this:scalar(\$q1)*this:scalar(\$q2) + v:dot(this:vector(\$q1), this:vector(\$q2))
};

(:~
: cross()
: Quaternion cross product
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return q1Xq2
:)
declare function this:cross(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
this:quaternion(0, this:vector(\$q1=>this:multiply(\$q2)))
};

(:~
: pow()
: Raise a quaternion to the given power.
:
: @param \$q: the quaternion
: @param \$pow: the power to raise it to
: @return q^pow
:)
declare function this:pow(
\$q as xs:double*,
\$pow as xs:double
) as xs:double*
{
if (this:is-scalar(\$q) and (this:ps(\$q) >= 0 or \$pow >= 1)) then (
this:as-quaternion(math:pow(this:ps(\$q), \$pow))
) else
switch (\$pow)
case 0 return \$this:one
case 1 return \$q
case 2 return
this:quaternion(
this:ps(\$q)*this:ps(\$q) - this:pi(\$q)*this:pi(\$q) - this:pj(\$q)*this:pj(\$q) - this:pk(\$q)*this:pk(\$q),
2*this:ps(\$q)*this:pi(\$q),
2*this:ps(\$q)*this:pj(\$q),
2*this:ps(\$q)*this:pk(\$q)
)
case 3 return
let \$a := this:ps(\$q)
let \$b := this:pi(\$q)
let \$c := this:pj(\$q)
let \$d := this:pk(\$q)
return (
this:quaternion(
math:pow(\$a,3) - 3*\$b*\$b*\$a - 3*\$c*\$c*\$a - 3*\$d*\$d*\$a,
3*\$a*\$a*\$b - math:pow(\$b,3) - 3*\$c*\$c*\$a - 3*\$d*\$d*\$a,
3*\$a*\$a*\$c - math:pow(\$c,3) - 3*\$b*\$b*\$c - 3*\$d*\$d*\$c,
3*\$a*\$a*\$d - math:pow(\$d,3) - 3*\$b*\$b*\$d - 3*\$c*\$c*\$d
)
)
default return (
this:exp(this:log(\$q)=>v:times(\$pow))
)
};

(:~
: log()
: Natural log of quaternion.
:
: @param \$q: the quaternion
: @return log(q)
:)
declare function this:log(
\$q as xs:double*
) as xs:double*
{
let \$m := this:magnitude(\$q)
return (
this:quaternion(
math:log(\$m),
(this:vector(\$q)=>v:normalize())=>
v:times(math:acos(this:scalar(\$q) div \$m))
)
)
};

(:~
: exp()
: Exponential of q
:
: @param \$q: the quaternion
: @return e^q
:)
declare function this:exp(
\$q as xs:double*
) as xs:double*
{
let \$ea := math:exp(this:scalar(\$q))
let \$nv := v:magnitude(this:vector(\$q))
return (
this:quaternion(
\$ea * math:cos(\$nv),
(this:vector(\$q)=>v:normalize())=>v:times(\$ea * math:sin(\$nv))
)
)
};

(:~
: sqrt()
: Principal square root of quaternion.
:
: @param \$q: the quaternion
: @return √q
:)
declare function this:sqrt(
\$q as xs:double*
) as xs:double*
{
if (this:is-scalar(\$q)) then (
if (this:ps(\$q) >= 0) then this:as-quaternion(math:sqrt(this:ps(\$q)))
else this:quaternion(0, math:sqrt(abs(this:ps(\$q))), 0, 0)
) else (
this:quaternion(
math:sqrt((this:magnitude(\$q) + this:ps(\$q)) div 2),
v:normalize(this:vector(\$q))=>
v:times(math:sqrt((this:magnitude(\$q) - this:ps(\$q)) div 2))
)
)
};

(:~
: sin()
: Sine of quaternion. Warning: some tests involving q:sin() are still failing
: (although many are not). Still working it out.
:
: @param \$q: the quaternion
: @return sin(q)
:)
declare function this:sin(
\$q as xs:double*
) as xs:double*
{
(:
: sin(q) = -1/2 n(e^(q n) - e^(-q n)) |v(q)|!=0; sin(s) otherwise
: where n = v(q)/|v(q)| i.e. normalized vector(q)
: which is -n sinh(qn)
:)
if (abs(v:magnitude(this:vector(\$q))) > \$config:ε) then (
let \$nv := this:quaternion(0, this:vector(\$q)=>v:normalize())
return this:minus(\$nv)=>this:multiply( this:sinh(\$q=>this:multiply(\$nv)) )
) else (
this:as-quaternion(math:sin(this:scalar(\$q)))
)
};

(:~
: cos()
: Cosine of quaternion.
:
: @param \$q: the quaternion
: @return cos(q)
:)
declare function this:cos(
\$q as xs:double*
) as xs:double*
{
(:
: cos(q) = 1/2 (e^(q n) + e^(-q n)) |q|!=0; cos(s) otherwise
: where n = v(q)/|v(q)| i.e. normalized vector(q)
: which is cosh(qn)
:)
if (abs(v:magnitude(this:vector(\$q))) > \$config:ε) then (
let \$nv := this:quaternion(0, this:vector(\$q)=>v:normalize())
return this:cosh( \$q=>this:multiply(\$nv) )
) else (
this:as-quaternion(math:cos(this:scalar(\$q)))
)
};

(:~
: sinh()
: Hyperbolic sine.
:
: @param \$q: the quaternion
: @return sinh(q)
:)
declare function this:sinh(\$x as xs:double*) as xs:double*
{
(this:exp(\$x)=>this:sub( this:exp(this:minus(\$x)) ))=>this:times(1 div 2.0)
};

(:~
: cosh()
: Hyperbolic cosine.
:
: @param \$q: the quaternion
: @return cosh(q)
:)
declare function this:cosh(\$x as xs:double*) as xs:double*
{
};

(:~
: times()
: Multiply a quaternion by a constant.
:
: @param \$q: the quaternion
: @param \$k: the constant to multiply by
: @return k*q
:)
declare function this:times(
\$q as xs:double*,
\$k as xs:double
) as xs:double*
{
\$q=>v:times(\$k)
};

(:~
: plus()
: Add a constant to a quaternion. (i.e. add quaternion (\$k, 0, 0, 0))
:
: @param \$q: the quaternion
: @param \$k: the constant to add
: @return q+k
:)
declare function this:plus(
\$q as xs:double*,
\$k as xs:double
) as xs:double*
{
};

(:~
: divide()
: Divide two quaternions q1/q2 as q1*q2^-1
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return q1/q2
:)
declare function this:divide(
\$q1 as xs:double*,
\$q2 as xs:double*
) as xs:double*
{
\$q1=>this:multiply(this:reciprocal(\$q2))
};

(:~
: magnitude()
: Magnitude of a quaternion; aka modulus()
:
: @param \$q: the quaternions
: @return length of quaternion vector
:)
declare function this:magnitude(
\$q as xs:double*
) as xs:double
{
v:magnitude(\$q)
};

(:~
: modulus()
: Same as magnitude(), name you see in math papers
:
: @param \$q: the quaternions
: @return modulus
:)
declare function this:modulus(
\$q as xs:double*
) as xs:double
{
v:magnitude(\$q)
};

(:~
: normalize()
: Normalize the quaternion so magnitude is 1. aka sgn()
: Fails for (0,0,0,0)
:
: @param \$q: the quaternion to normalize
: @returns normalized quaternion
: @error Divide by zero for the origin
:)
declare function this:normalize(\$q as xs:double*) as xs:double*
{
let \$l := this:magnitude(\$q)
return (
this:quaternion(
this:ps(\$q) div \$l,
this:pi(\$q) div \$l,
this:pj(\$q) div \$l,
this:pk(\$q) div \$l
)
)
};

(:~
: sgn()
: Same as normalize(); name you see in math papers
: Doesn't speak to me.
:
: @param \$q: the quaternion to normalize
: @returns normalized quaternion
: @error Divide by zero for the origin
:)
declare function this:sgn(\$q as xs:double*) as xs:double*
{
this:normalize(\$q)
};

(:~
: conjugate()
: Return the quaternion conjugate
:
: @param \$q: quaternion
: @return conj(q)
:)
declare function this:conjugate(\$q as xs:double*) as xs:double*
{
this:quaternion(this:ps(\$q), -this:pi(\$q), -this:pj(\$q), -this:pk(\$q))
};

(:~
: reciprocal()
: The reciprocal of q: q^-1 = conj(q)/||q||²
:
: @param \$q: quaternion
: @return q^-1
:)
declare function this:reciprocal(\$q as xs:double*) as xs:double*
{
(: Avoid the extra sqrt :)
let \$n2 :=
this:ps(\$q)*this:ps(\$q) +
this:pi(\$q)*this:pi(\$q) +
this:pj(\$q)*this:pj(\$q) +
this:pk(\$q)*this:pk(\$q)
return this:conjugate(\$q)=>this:times(1 div \$n2)
};

(:~
: abs()
: The absolute value of z
:
: @param \$q: quaternion
: @return abs(q)
:)
declare function this:abs(\$q as xs:double*) as xs:double*
{
this:quaternion(\$q=>v:abs())
};

(:~
: quote()
: Return a string value of the quaternion
: @param \$q: quaternion
: @return string
:)
declare function this:quote(\$qs as xs:double*) as xs:string
{
string-join(
let \$chunks :=
for \$q at \$i in \$qs return (
if (\$i mod 4 = 0) then (
array {\$qs[\$i - 3], \$qs[\$i - 2], \$qs[\$i - 1], \$q}
) else ()
)
for \$chunk in \$chunks
let \$q := \$chunk?*
return (
let \$s := this:scalar(\$q)
let \$v := this:vector(\$q)
let \$ps := if (\$s != 0 or (every \$c in \$v satisfies \$c = 0)) then string(\$s) else ()
let \$pi :=
if (\$v[1] = 0) then ()
else if (\$v[1] = -1) then "-i"
else if (empty(\$ps)) then (if (\$v[1]=1) then "i" else \$v[1]||"i")
else if (\$v[1] < 0) then \$v[1]||"i"
else if (\$v[1] = 1) then "+i"
else "+"||\$v[1]||"i"
let \$pj :=
if (\$v[2] = 0) then ()
else if (\$v[2] = -1) then "-j"
else if (empty((\$ps, \$pi))) then (if (\$v[2]=1) then "j" else \$v[2]||"j")
else if (\$v[2] < 0) then \$v[2]||"j"
else if (\$v[2] = 1) then "+j"
else "+"||\$v[2]||"j"
let \$pk :=
if (\$v[3] = 0) then ()
else if (\$v[3] = -1) then "-k"
else if (empty((\$ps, \$pi, \$pj))) then (if (\$v[3]=1) then "k" else \$v[3]||"k")
else if (\$v[3] < 0) then \$v[3]||"k"
else if (\$v[3] = 1) then "+k"
else "+"||\$v[3]||"k"
return (
if (\$s=0 and exists((\$pi,\$pj,\$pk)))
then \$pi||\$pj||\$pk
else \$ps||\$pi||\$pj||\$pk
)
)
,
" "
)
};

(:~
: unquote()
: Parse a string value of a quaternion to produce that quaternion.
: Number is of the form 33+12i+6j+8k
:
: @param \$qs: the string representation of the number
: @return quaternion
:)
declare function this:unquote(\$qs as xs:string) as xs:double*
{
let \$analysis :=
analyze-string(
replace(\$qs, " ", ""),
"([+-]?((\d+([.]\d)?[ijk]?)|(\d*[.]\d+[ijk?])|([ijk])))"
)
let \$parts := \$analysis/fn:match
return (
if (count(\$parts) > 4 or exists(\$analysis/fn:non-match))
else (
let \$six :=
(for \$part at \$i in \$parts where not(matches(\$part, "[ijk]")) return \$i)
let \$iix :=
for \$part at \$i in \$parts where ends-with(\$part, "i") return \$i
let \$jix :=
for \$part at \$i in \$parts where ends-with(\$part, "j") return \$i
let \$kix :=
for \$part at \$i in \$parts where ends-with(\$part, "k") return \$i
let \$number :=
function(\$ix as xs:integer?, \$v as xs:string) as xs:double {
if (empty(\$ix)) then 0 else (
let \$s := string(\$parts[\$ix])
let \$prefix := substring-before(\$s, \$v)
return (
if (empty(\$s)) then 0
else if (\$prefix=("","+")) then 1
else if (\$prefix="-") then -1
else number(\$prefix)
)
)
}
return (
if (count(\$six) > 1 or count(\$iix) > 1 or count(\$jix) > 1 or count(\$kix) > 1)
else (
this:quaternion(
if (exists(\$six)) then number(string(\$parts[\$six])) else 0,
\$number(\$iix, "i"),
\$number(\$jix, "j"),
\$number(\$kix, "k")
)
)
)
)
)
};

(:~
: same()
: Are the two quaternions the same?
:
: @param \$q1: one quaternion
: @param \$q2: other quaternion
: @return true() if the values are the same
:)
declare function this:same(\$q1 as xs:double*, \$q2 as xs:double*) as xs:boolean
{
this:ps(\$q1)=this:ps(\$q2) and
this:pi(\$q1)=this:pi(\$q2) and
this:pj(\$q1)=this:pj(\$q2) and
this:pk(\$q1)=this:pk(\$q2)
};

(:~
: decimal()
: Cast the values to a values with the given number of digits.
:
: @param \$qs: the quaternions
: @param \$digits: how many digits to preserve
: @return complex numbers
:)
declare function this:decimal(\$qs as xs:double*, \$digits as xs:integer) as xs:double*
{
v:decimal(\$qs, \$digits)
};```