---
title: "Matrices with quaternion or octonion entries: class `onionmat` in the `onion` package"
author: "Robin K. S. Hankin"
output: html_vignette
bibliography: onionpaper.bib
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{onionmat}
  %\usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
options(rmarkdown.html_vignette.check_title = FALSE)
library("onion")
set.seed(1)
```

```{r out.width='20%', out.extra='style="float:right; padding:10px"',echo=FALSE}
knitr::include_graphics(system.file("help/figures/onion.png", package = "onion"))
```

## Introduction

To cite the `onion` package in publications please use
@hankin2006_onion.  This short document shows how matrices with
octonion or quaternion entries may be created and manipulated with the
`onionmat` class of the `onion` R package.

## `onionmat` objects

A good place to start is function `romat()`, which creates a simple
`onionmat` object:

```{r}
set.seed(0)
options(digits=2)
romat()
```

This illustrates many features of the package.  An `onionmat` object
has two slots.  The first slot, `x`, is an onion [a vector of
quaternions or octonions] and the second, `M`, an integer matrix
which is used to store attributes such as dimensions and dimnames.
The elements of `M` and `d` are in bijective correspondence; thus
element `[b,AR]` is number 17 and this is seen to be approximately
$-0.81 + 0.24i -1.43j +0.37k$.  Most R idiom will work with such
objects, here is a brief sample.

```{r}
A <- matrix(rquat(21),7,3)  # matrix() calls onion::onionmat()
A
```

See above how object `A` has no rownames or colnames and the defaults are used.  We may extract components:

```{r}
A[1,]
```

above, the resulting object is an onion but we may retain the onionmat character using `drop`:

```{r}
A[1,,drop=FALSE]
```


The extraction methods operate as expected:

```{r}
Re(A)
k(A)
```


(above, the matrices returned are numeric).  Also replacement methods work as expected:

```{r}
j(A) <- -1
A
```

Some of the summary methods work:

```{r}
sum(A)
```

## Matrix multiplication

Matrix multiplication is implemented.


```{r}
A <- matrix(rquat(21),3,7)
umbral <- state.abb[1:7]
rownames(A) <- letters[1:3]
colnames(A) <- umbral

B <- matrix(rquat(28),7,4)
rownames(B) <- umbral
colnames(B) <- c("H","He","Li","Be")

A %*% B
```


However, it is often preferable (but no faster in this case) to use
functions such as `cprod()` and `tcprod()`:


```{r}
C <- matrix(rquat(14),7,2)
rownames(C) <- umbral
colnames(C) <- month.abb[1:2]
cprod(B,C)
```

and indeed the single-argument versions work as expected:

```{r}
tcprod(A) - A %*% ht(A)
```


## Octonions

Here I consider $3\times 3$ Hermitian octonionic matrices which are
known be a special case of a Jordan algebra (see the `jordan` package
for more systematic Jordan algebra R functionality).

```{r}
x <- cprod(matrix(roct(12),4,3))
x
```

We see that `x` is Hermitian symmetric:

```{r}
max(Mod(Im(x+t(x))))
```

[that is, the imaginary components of symmetrically placed elements
are mutually conjugate].  We may verify that $3\times 3$ matrices form
a Jordan algebra under the composition rule $A\circ B=(AB+BA)/2$
[juxtaposition indicating regular matrix multiplication]; the identity is 

\[(xy)(xx) = x(y(xx)).\]

First we define the Jordan product:

```{r}
`%o%` <- function(x,y){(x%*%y + y%*%x)/2}
```

then create a couple of random Hermitian octonionic matrices:

```{r}
x <- cprod(matrix(roct(12),4,3))
y <- cprod(matrix(roct(12),4,3)) # x and y are 3-by-3 Hermitian octonionic matrices
```

We first verify numerically that the Jordan product of two Hermitian
symmetric matrices is Hermitian:

```{r}
jj <- x %o% y
max(Mod(Im(jj+t(jj))))
```

then verify the Jordan identity:

```{r}
LHS <- (x %o% y) %o% (x %o% x)
RHS <- x %o% (y %o% (x %o% x))
max(Mod(LHS-RHS))  # zero to numerical precision
```

showing that the Jordan identity is satisfied, up to a small numerical
tolerance.  However, $4\times 4$ octonionic matrices do not satisfy
the Jordan identity:

```{r}
x <- cprod(matrix(roct(16),4,4))
y <- cprod(matrix(roct(16),4,4)) # x and y are 4-by-4 Hermitian octonionic matrices
LHS <- (x %o% y) %o% (x %o% x)
RHS <- x %o% (y %o% (x %o% x))
max(Mod(LHS-RHS))  # miles off
```


## Note on the print method

Above we have been using the default print method but it is possible
to use a more compact notation, which is useful if a matrix has many
short entries.


```{r}
o <- rsomat()
o
options(show_onionmats_in_place = TRUE)
o
```