Introduction

The following library aims to deliver a simple method to construct hypercomplex numbers from any of the Cayley-Dickson algebras. It supports calculations in an arbitrary-precise arithmetic as well as encryption/decryption procedures for public-key lattice-based cryptosystems in high-dimensional algebras of truncated polynomial rings. The library is dedicated mostly to computational mathematicians and computational scientists whose focus is (post-quantum) cryptography and precise computation. As a header-only C++ template code it's greatest advantage is the combination of speed, generic programming and convenience for the end user, making it well suited for wide range of computationally challenging projects.

Installation

The following is a header-only library, meaning that the easiest way to use it is to copy the core hpp files alongside your main program and include it into the code with the directive:

#include "Hypercomplex.hpp"

Remember to specify a proper langauge standard and optimisation level for the compiler, as in the command below:

g++ -O2 --std=c++17 main.cpp -o main

For computations based on high precision arithmetics it is essential to install the MPFR library first. This should be rather straightforward with the following commands:

# Linux (apt-get package manager)
sudo apt-get update -y
sudo apt-get install -y libmpfr-dev
 
# macOS (Homebrew package manager)
brew update
brew install mpfr

It is then essential to provide a preprocessor macro specifying the use of this library at the compile time as well as linking with GNU MP and GNU MPFR as in the command below:

g++ -O2 -DUSEMPFR=1 --std=c++17 main.cpp -o main -lmpfr -lgmp

Alternatively, if you work with conda environments this library is installable with:

conda install -c angrymaciek hypercomplex

However, please note that during compilation time it might be necessary to provide additional flags which point to the headers and other libs for the linker.

The command above would be adjusted to:

g++ -O2 -DUSEMPFR=1 --std=c++17 -I$CONDA_PREFIX/include -L$CONDA_PREFIX/lib main.cpp -o main -lmpfr -lgmp

Brief overview

The following section demonstrates general functionality and behaviour of the library. For the full unit test suite please refer to this file.

Please note that throughout this whole documentation many links may point to the main template Hypercomplex class instead of its specialisations. This is because the doxygen engine unfortunately cannot distinguish between template classes of the same name properly. Always check the description carefully while accessing links at this page.

Let us construct two hypercomplex numbers from an algebra obtained with the Cayley-Dickson process. For simplicity of the presentation we will focus on quaternions.

\( H_1=(1,0,-0.5,5)\)

\( H_2=(-2,-4,-6,0)\)

double arr1[4] = {1.0,0.0,-0.5,5.0};
Hypercomplex<double, 4> H1(arr1);
std::cout << "H1 = " << H1 << std::endl;
 
double arr2[4] = {-2.0,-4.0,-6.0,0.0};
Hypercomplex<double, 4> H2(arr2);
std::cout << "H2 = " << H2 << std::endl;

The code above results in:

H1 = 1 0 -0.5 5
H2 = -2 -4 -6 0

For every hypercomplex number we may extract its' real as well as imaginary part:

\(Re(H) := (H^{(0)}, 0, 0, 0, \dotsc)\)

\(Im(H) := (0, H^{(1)}, H^{(2)}, H^{(3)}, \dotsc)\)

Therefore the commands:

std::cout << "Re(H1) = " << Re(H1) << std::endl;
std::cout << "Im(H1) = " << Im(H1) << std::endl;

yield:

Re(H1) = 1 0 0 0
Im(H1) = 0 0 -0.5 5

In case you already forgot what the dimensionality of our objects is, don't worry - we got your back. Asking the right questions...

std::cout << "dim(H1) = " << H1._() << std::endl;

...leads to right answers:

dim(H1) = 4

It would be very nice if we could access the components directly... But we can!

std::cout << "H2(2) = " << H2[2] << std::endl;

results in:

H2(2) = -6

Do you wish you could represent your objects in different systems? Nothing easier than that. Embedding lower-dimensional hypercomplex numbers into higher dimensional algebras is a piece of cake!

std::cout << "Oct(H1) = " << H1.expand<8>() << std::endl;

The template method above creates a new class instance, printed below:

Oct(H1) = 1 0 -0.5 5 0 0 0 0

For every hypercomplex number we might calculate its' Euclidean norm:

\(||H||_2 := \sqrt{H^{(0)}\times H^{(0)} + \dotsc + H^{(n-1)} \times H^{(n-1)} }\)

Luckily, there is a function dedicated to this operation:

std::cout << "||H2|| = " << H2.norm() << std::endl;

It returns value of the same type as the template type for the base class:

||H2|| = 7.48331

Having defined a norm for every non-zero hypercomplex number we may calculate its' inverse according to the following formula: \(H^{-1} = \frac{\bar{H}}{||H||_2^2}\)

Calling a class method...

std::cout << "H2^-1 = " << H2.inv() << std::endl;

...creates a new class instance:

H2^-1 = -0.0357143 0.0714286 0.107143 -0

Arithmetic operations on hypercomplex numbers are not complicated at all.

• Addition is carried out over respective elements: \(\forall_{H_A, H_B}: [H_A+H_B]^{(i)} := H_A^{(i)} + H_B^{(i)}\)
• Subtraction is carried out over respective elements: \(\forall_{H_A, H_B}: [H_A-H_B]^{(i)} := H_A^{(i)} - H_B^{(i)}\)
• A general formula for multiplication of hypercomplex numbers may be derived as follows: Let \(H_A\) and \(H_B\) be elements from a Cayley-Dickson algebra of dimension \(2^n\). Both numbers may be interpreted as ordered pairs of elements from a \(2^{(n-1)}\)-dimensional algebra: \(H_A = (a,b)\) and \(H_B = (c,d)\). Such a representation yields a recursive multiplication algorithm: \(H_A \times H_B = (a,b)(c,d) := (ac-\bar{d}b,da+b\bar{c})\). (Multiplication of hypercomplex numbers is indeed implemented as a recursive operator with the base condition of multiplying real numbers.) Disclaimer: Various distinct definitions of the multiplication formula exist: here, here or here. They all lead to a proper norm upon multiplying a number with it's conjugate and the choice is arbitrary (feel free to modify the code, if needed). In addition to the overloaded operator, hypercomplex multiplication is also available as a class method. Depending on the dimenion of algebra and compiler's optimisation settings it could be 5-50x faster than the overloaded operator.
• Knowing that inverse elements of hypercomplex numbers exist a division operation is implementat as a multiplication with an inverse of the right operand: Let \(H_A\) and \(H_B\) be elements from a Cayley-Dickson algebra of dimension \(2^n\), \(H_B \neq 0\). \(\frac{H_A}{H_B} := H_A \times H_B^{-1}\) Notice the order of the operands, as commutativity is no longer a given.

To test these operations we may execute the code below:

std::cout << "H1 + H2 = " << H1 + H2 << std::endl;
std::cout << "H1 - H2 = " << H1 - H2 << std::endl;
std::cout << "H1 * H2 = " << H1 * H2 << std::endl;
std::cout << "H1 / H2 = " << H1 / H2 << std::endl;

Obtained output should match:

H1 + H2 = -1 -4 -6.5 5
H1 - H2 = 3 4 5.5 5
H1 * H2 = -5 26 -25 -12
H1 / H2 = 0.0178571 -0.464286 0.482143 -0.142857

Moreover, one can easily raise a hypercomplex number to a natural power: \((H^n, n\in N)\). This operation is implemented simply as an iterative multiplication:

\(H^n := (\dotsc((H \times H) \times H)\dotsc)\)

Please note that the products are evaluated from LHS to RHS. However, all Cayley-Dickson algebras are power-associative therefore it does not really matter.

Calling:

std::cout << "H2^4 = " << (H2^4) << std::endl;

Produces:

H2^4 = 1472 -1536 -2304 0

Last, but not least: our little cherry on top. A very special function linked to the magical Euler number - hypercomplex exponentiation:

\(e^H = e^{Re(H)+Im(H)} := e^{Re(H)} \times (cos||Im(H)||_2 + \frac{Im(H)}{||Im(H)||_2} \times sin||Im(H)||_2)\)

Regardless of the algebra hypercomplex numbers in the equation above are multiplied by scalars, therefore associativity and commutativity still holds for these formulas. For that reason the exponentiation function is highly optimized, implemented efficiently with as few operations and variables as possible.

Finally, we validate it with:

std::cout << "e^H1 = " << exp(H1) << std::endl;

and get:

e^H1 = 0.83583 -0 0.257375 -2.57375

Arbitrary-precision arithmetic

Calculations on MPFR types are availabla via partial template specialisation (meaning that the hpp file already includes the necessary header and the user should not include it again). It is strongly advised to read the library manual beforehand (sections 1 and 4 as a must-read). Please do not mix this library with another MPFR-dependant code in the same translation unit.

We start with setting a global precision for all objects.

set_mpfr_precision(200);

What follows is an initialization of 8 MPFR variables in a single array:

mpfr_t A[8];
mpfr_init2(A[0], get_mpfr_precision());
mpfr_init2(A[1], get_mpfr_precision());
mpfr_init2(A[2], get_mpfr_precision());
mpfr_init2(A[3], get_mpfr_precision());
mpfr_init2(A[4], get_mpfr_precision());
mpfr_init2(A[5], get_mpfr_precision());
mpfr_init2(A[6], get_mpfr_precision());
mpfr_init2(A[7], get_mpfr_precision());
mpfr_set_d(A[0], 1.5, MPFR_RNDN);
mpfr_set_d(A[1], 2.5, MPFR_RNDN);
mpfr_set_d(A[2], 0.0, MPFR_RNDN);
mpfr_set_d(A[3], -1.5, MPFR_RNDN);
mpfr_set_d(A[4], 0.5, MPFR_RNDN);
mpfr_set_d(A[5], -0.5, MPFR_RNDN);
mpfr_set_d(A[6], -0.5, MPFR_RNDN);
mpfr_set_d(A[7], -1.5, MPFR_RNDN);

Let us create an octonion composed of these numbers:

Hypercomplex<mpfr_t, 8> Hx(A);

To print the first element of our number after raising it to the 30-th power we need to use library specific function:

std::cout << "Hx^30 = ";
mpfr_out_str(stdout, 10, 0, (Hx^30)[0], MPFR_RNDN);
std::cout << std::endl;

In result we obtain:

Hx^30 = -1.1841044160704622190868522692471742630004882812500000000000000e17

After all the calculations it is essential to clear constructed objects from the memory manually:

mpfr_clear(A[0]);
mpfr_clear(A[1]);
mpfr_clear(A[2]);
mpfr_clear(A[3]);
mpfr_clear(A[4]);
mpfr_clear(A[5]);
mpfr_clear(A[6]);
mpfr_clear(A[7]);

Also, following that, to call a wrapper function which cleans all internally-reserved memory:

clear_mpfr_memory();

All the code specified up to this point may be executed upon compilation of this source code.

Cryptographic Application

The main feature of Hypercomplex is an implementation of cryptographic operations as in NTRU cryptosystem but generalized on an arbitrary-high-dimensional algebras generated with the Cayley-Dickson construction. The library provides additional helper class for truncated polynomials (with necessary operators e.g. modular reduction, convolution-multiplication) and a partial template specialisation for hypercomplex numbers based on these objects.

Briefly, let \(N, p, q\) denote three selected primes such that \(p << q\). The former will become a multiplicative order whereas the two latter will mark characteristics of the structures we will work on. Let \(R = \mathbb{Z}[x] / (x^N - 1)\) be the underlying polynomial quotient ring ( \(R_p, R_q\) are modular structures, modulo \(p,q\) respectively) and \(D\) mark the dimension of the algebra we operate in: \(A = \{x_0 + \sum^{D-1}_{i=1} x_i \cdot e_i | x_i \in R\} \), where all \(e_i\) are imaginary basis elements (similarly for modular algebras: \(A_p, A_q\)).

Mathematical derivations for these structures and their operations are analogous to those presented for QTRU and OTRU.

Imagine that Alice wishes to send a message to Bob. Bob chooses \(F \in A\) such that \(F_p^{-1} \in A_p\) and \(F_q^{-1} \in A_q\) exist, then generates his public key for the cryptosystem:

\( H = F_q^{-1} \times G \mod q\)

Alice chooses a "blinding element" \(\Phi \in A_q\) and encrypts her top secret message \(M\) as:

\( E = (p \cdot H \times \Phi + M) \mod q\)

Which Bob then decrypts with the following operations:

\( C_1 = ((F \times E) \times F) \mod q\)

\( C_2 = C_1 \mod p\)

\( C_3 = (F_p^{-1} \times (C_2 \times F_p^{-1})) \mod p\)

If the decryption was successfull Bob receives \( C_3 = M\) (up to coefficients' centered lift)

Please remember that lattice-based cryptography is always burdened with a chance of decryption failure due to incorrect recovery of polynomial's coefficients at the centered lift step.

With a simple example of how to create a hypercomplex number based on truncated polynomials...

int64_t X_array1[] = {2, 1, 1, 0, 2};
int64_t X_array2[] = {2, 2, 2, 1, 2};
Polynomial<5> X_polynomial1(X_array1);
Polynomial<5> X_polynomial2(X_array2);
Polynomial<5> X_coefficients[] = { X_polynomial1, X_polynomial2 };
Hypercomplex<Polynomial<5>, 2> X(X_coefficients);

... the scheme presented above may be implemented with the following functions of the library:

# Having defined: p, q, MaxDeg, dim, F, G, M, PHI
 
# PUBLIC KEY
CenteredLift(&F, p);
CenteredLift(&G, p);
Hypercomplex<Polynomial<MaxDeg>, dim> H = PUBLICKEY(F, G, q);
 
# ENCRYPTION
CenteredLift(&PHI, p);
Hypercomplex<Polynomial<MaxDeg>, dim> E = ENCRYPT(H, M, PHI, p, q);
 
# DECRYPTION
Hypercomplex<Polynomial<MaxDeg>, dim> C = DECRYPT(F, E, p, q);

Remarkably, for a cryptosystem based on an algebra of \(D \geq 16\) dimensions \(F\) needs to contain at most one \( x_i \in R | x_i \neq 0\). This is because sedonions (and higher) are not associative, thus the decryption process will only be possible for a specific, reduced subset of private keys.

Cryptographic applications of Hypercomplex have been extensively tested in the test case: Cryptosystem based on Cayley-Dickson Algebras of the following file.

All tests underlying Figure 1 in the publication are available here.