Thanks to Gregg Stiesberg for bug report.

§1a. This report offers a simple C++14 implementation of the multicomplex numbers. Everyone may use, modify, and republish the source code.

Loosely speaking, the multicomplex numbers (abbreviated MC) are complex numbers whose elements are themselves complex numbers, and this recursion can be carried to arbitrary, but finite, depth. For a formal mathematical development, start with these:

• Let i₁, i₂, i₃ … be a sequence of numbers such that i_n² = −1.
• For the operations of addition and multiplication, postulate that all the i_n's and real numbers are commutative and associative among one another, and that multiplication is distributive over addition. However, neither the sum nor product of two or more distinct i_n's can be simplified.

Now define symbols ℂ_n:

• ℂ₀ is the well-known field of the real numbers.
• ℂ₁ is the well-known field of the complex numbers over the reals. Each number in ℂ₁ is an ordered pair (x, y) of numbers taken from ℂ₀, and is interpreted as x + yi₁.
• Recursively, ℂ_n+1 is the ring (but not field) of the multicomplex numbers over ℂ_n. Each number in ℂ_n+1 is an ordered pair (x, y) of numbers taken from ℂ_n, and is interpreted as x + yi_n+1.

The subscript of ℂ is the order of the multicomplex number.

Often, x is termed the real part and y the imaginary part. With the recursion, care is advised in using these terms because the real part typically contains imaginary parts as subelements, and vice versa.

§1b. Multicomplex arithmetic is based on that of ordinary complex numbers. In the following examples two MC's will be employed: w, z ∈ ℂ_n+1. These can respectively be written w = u + vi_n+1 and z = x + yi_n+1; where u, v, x, y ∈ ℂ_n.

Being of the same order, w = z if and only if u = x and v = y. With two numbers of different order, equality can sometimes be achieved by zero-padding in the imaginary part. Given s ∈ ℂ_n, construct t ∈ ℂ_n+1 such that t = s + 0i_n+1. (This can be repeated to obtain higher orders.) Then t is deemed equal to s. For a commonplace example of zero-padding, recall that a ℂ₀ number is often regarded as a ℂ₁ number whose imaginary part is zero.

Some familiar operations are defined as follows:

• w + z = (u + x) + (v + y) i_n+1
• w − z = (u − x) + (v − y) i_n+1
• wz = (ux − vy) + (uy + vx) i_n+1

For r a real number, wr equals ur + vri_n+1 because r can be zero-padded into an MC of appropriate order, and then the above multiplication rule can be applied. More broadly, had w and z been of any two different orders, the number of smaller order could have been zero-padded repeatedly until the two numbers had equal order. Then addition, subtraction, or multiplication could take place.

The conjugate of an MC, here written with a prefix tilde, negates the second part. In symbols, ~w = u − vi_n. Note that all the real and imaginary elements and subelements within v are negated.

The norm of an MC, which is always a real number, is defined recursively:

• If z ∈ ℂ₀, then norm (z) = z²
• Otherwise, norm (z) = norm (x) + norm (y)

The folding operation multiplies a number by its conjugate. With w ∈ ℂ_n+1, fold (w) = w(~w) = (~w)w ∈ ℂ_n. Meanwhile, the fold of a ℂ₀ number is simply its square. Sometimes the fold of a nonzero number turns out to be zero; for example, (1 + i₁i₂)(1 − i₁i₂) = 0. This illustrates that the norm is not always preserved under multiplication, thus the MC's are not a composition algebra. After sufficiently many foldings, an MC becomes a real number; if it turns out to be zero, the MC is abient (from a Latin participle meaning "going away", and not to be confused with ambient).

There is a recursive algorithm that will produce the quotient in division when the divisor is not abient:

1. To divide an MC by a real number (i.e. ℂ₀), multiply the MC by the reciprocal of the real. Then the algorithm is complete.
2. To divide an MC (the dividend) by another MC (the divisor) when the divisor's order is greater than zero:
   • Multiply the dividend by the conjugate of the divisor to obtain the new dividend.
   • Multiply the divisor by the conjugate of the divisor to obtain the new divisor. This effects a folding in the divisor, so that the order of the new divisor is one less than the order of the old divisor.
3. Repeat step two until the new divisor is a real number. Then go to step one.

§1c. In higher-order multicomplex numbers, products like i₂i₄i₅ can appear, and can become notationally awkward. An alternative employs the letter j and a pattern that resembles binary numbers:

j1 = i1		j5 = i3 i1
j2 = i2		j6 = i3 i2
j3 = i2 i1		j7 = i3 i2 i1
j4 = i3		j8 = i4

The subscript of j is calculated by addition:

in	addend
i1	1
i2	2
i3	4
in	2n−1

The symbol j₀ can be identified with the real number unity, and might be used for emphasis when none of the i_n apply.

§2a. The C++ code supplied with this page uses a custom-built template class for the multicomplex numbers. It allows many different numerical types for use as the real numbers, although double-precision floating point is the usual choice for many C++ programmers.

Here are some highlights of class multicomplex:

template <typename elem, int order>
class multicomplex {
public:
    using one_less = multicomplex<elem, order−1>;

    one_less real, imag; // these are the only data

    multicomplex (
        one_less const & r,
        one_less const & i
    );
    static multicomplex random ();

    static multicomplex const zero;
    static multicomplex const real_unit;
    static multicomplex const imag_unit;

    multicomplex operator+ () const;
    multicomplex operator~ () const;
    multicomplex operator− () const;

    template <int o_order>
    operator+ (multicomplex<elem, o_order> const & o) const;

    template <int o_order>
    operator− (multicomplex<elem, o_order> const & o) const;

    template <int o_order>
    operator* (multicomplex<elem, o_order> const & o) const;

    multicomplex operator* (elem const & o) const;

    template <int o_order>
    multicomplex operator/ (multicomplex<elem,o_order> const & o) const;

    multicomplex operator/ (elem const & o) const;

    one_less fold () const;

    template <int x_order>
    multicomplex<elem,x_order> const & at_con (int const n) const;

    template <int x_order>
    multicomplex<elem,x_order> & at_var (int const n);

    elem norm () const;

    void view_i (ostream &) const;
    void view_j (ostream &) const;
    void view_k (ostream &) const;
};

Because each order of MC's is implemented as a separate C++ class type, the code uses numerous templates. Further, multicomplex<elem,0> is a specialization of multicomplex<elem,depth>, and is generally similar to it, handling the ultimate case.

The real and imag elements of multicomplex<elem,depth> are public because there was no special reason to make them private. Also of public visibility is the sole data element of multicomplex<elem,0>, which is named sole; it has the mathematical characteristics of a real number.

Function random is a convenient way to get some random-looking MC's for testing. It is enabled only when elem is double.

§2b. In writing a program, certain using declarations are helpful. Assuming double for the fundamental numerical type:

ℂ0	using mcd0 = multicomplex<double,0>;
ℂ1	using mcd1 = multicomplex<double,1>;
ℂ2	using mcd2 = multicomplex<double,2>;
ℂ3	using mcd3 = multicomplex<double,3>;
et cetera

One way to construct multicomplex objects is to use initializer_lists, often nested:

mcd0 const a   {+1};
mcd1 const b   {+1,-2};
mcd2 const c  {{+1,-2},{-3,+4}};
mcd3 const d {{{+1,-2},{-3,+4}},{{+5,-6},{+7,+8}}};

A lower-order MC might appear within the initializer_list of a higher-order MC, so that the following set of definitions is equivalent to the previous:

mcd0 const a {+1};
mcd1 const b {a,-2};
mcd2 const c {b,{-3,+4}};
mcd3 const d {c,{{+5,-6},{+7,+8}}};

There are three ways to display an MC:

• view_i writes the elements additively, using i_n symbols.
• view_j writes the elements additively, using j_n symbols.
• view_k writes the MC's in nested form, with minimal non-space characters.

view_i:
a:     + 1.00
b:    (+ 1.00 - 2.00*i1)
c:   ((+ 1.00 - 2.00*i1) + (- 3.00 + 4.00*i1)*i2)
d:  (((+ 1.00 - 2.00*i1) + (- 3.00 + 4.00*i1)*i2) + ((+ 5.00 - 6.00*i1) + (+ 7.00 + 8.00*i1)*i2)*i3)

view_j:
a: (+ 1.00*j0)
b: (+ 1.00*j0 - 2.00*j1)
c: (+ 1.00*j0 - 2.00*j1 - 3.00*j2 + 4.00*j3)
d: (+ 1.00*j0 - 2.00*j1 - 3.00*j2 + 4.00*j3 + 5.00*j4 - 6.00*j5 + 7.00*j6 + 8.00*j7)

view_k:
a:     1.00
b:    [1.00 -2.00]
c:   [[1.00 -2.00] [-3.00 4.00]]
d:  [[[1.00 -2.00] [-3.00 4.00]] [[5.00 -6.00] [7.00 8.00]]]

For output stream insertion, operator<< invokes view_k.

§2c. Functions named at_con and at_var retrieve a sub-MC from a MC. Either can be called on a constant MC, and will return a constant reference to a sub-MC. On the other hand, at_var can be called on only a variable (i.e. non-const) MC, and it will return a variable reference to a sub-MC.

This differs from the behavior of std::vector::at which comes in both constant and variable versions, because with a variable vector there is no straightforward way to call the constant version of at. A cautious programmer might want to do that in order to prevent accidentally modifying the vector.

The use of multicomplex::at is best explained with a lengthy example.

It would be possible to create random-access iterators in the style of the Standard Library.

§2d. Unary operator+, operator~, and operator- have unsurprising behavior.

Binary operator+, operator-, and operator* may be used with two operands of different order. The program gives the same result as if the lower-order operand were zero-padded to match the higher-order, but a more efficient method is used. Binary operator/ also works with two operands of different order, but zero-padding is irrelevant.

For quicker performance, binary operator* and operator/ also have versions where the second operand is an ordinary real number.

§2e. C++ does not automatically supply any conversion operations between two instantiations of the same template. This is true even if the two instantiations involve types between which C++ offers a built-in conversion. For instance:

    double d {5.6 };
    float  f {8.9f};
    d = f; // OKAY -- implicit built-in conversion

    multicomplex<double,2> dm {{5.1 , 5.2 },{5.3 , 5.4 }};
    multicomplex<float ,2> fm {{8.1f, 8.2f},{8.3f, 8.4f}};
    dm = fm; // ERROR -- a conversion function is needed

As an example of a function to convert the element type, the source code provides convert_elem.

§2f. Supplied are a few transcendental functions: circular and hyperbolic sine and cosine; and the exponential function. As these are single-valued, calculation is straightforward. Many of the familiar trigonometric identities are satisfied.

A future version of the program may address the inverses of these functions, which mathematically are more difficult to manage because they are multiple-valued, and further because in the complex case there is often no strong criterion for deciding which value should be desginated principal.

Even the square root function becomes extremely multiple-valued for higher orders. For instance, here are some square roots of unity:

ℂ0	ℂ1	ℂ2	ℂ3	ℂ4		ℂ5
±1	±1	±1 ±i2i1	±1 ±i2i1 ±i3i1 ±i3i2	±1 ±i2i1 ±i3i1 ±i3i2	±i4i1 ±i4i2 ±i4i3 ±i4i3i2i1	±1 ±i2i1 ±i3i1 ±i3i2	±i4i1 ±i4i2 ±i4i3 ±i5i1	±i5i2 ±i5i3 ±i5i4 ±i5i4i3i2	±i5i4i3i1 ±i5i4i2i1 ±i5i3i2i1 ±i4i3i2i1

§3a. Optional for the programmer is to use auxiliary class flat<elem,order>. A flat object contains the same information as a multicomplex object, but in flat the elements are stored in a std::vector rather than a recursive structure. Although the facilities of flat might have been housed in multicomplex, the latter is already complicated enough.

In flat are found the following operations which lose no information and hence are reversible:

• a constructor from multicomplex to flat
• a conversion operator from flat to multicomplex
• a conversion operator from flat to std::vector<elem> const &
• a constructor from std::vector<elem> const & to flat where vector must have the correct number of elements

A flat object and its corresponding multicomplex object must have the same elem type and order.

Also available are at functions allowing individual elements to be manipulated in much the same manner as std::vector. Too, flat has output functions view_i, view_j, and view_k mimicking those of multicomplex; and operator<< invokes view_k.

On the other hand, flat does not offer any direct way to merge two flats as the real and imaginary components of a higher-order MC, or to split an MC into its real and imaginary parts. Also omitted are arithmetic and the transcendental functions, which are more cleanly handled by way of multicomplex's recursive structure.

Class flat's main purpose is to serve as a conduit between multicomplex and the widely-used std::vector.

§4a. Because class multicomplex uses a recursive template structure, it can handle different orders of MC's, which need different amounts of memory, without performing dynamic memory allocation. This can reduce run time considerably.

MC's as supported by multicomplex must be balanced, in the sense that the real part of an MC must be of the same order as the corresponding imaginary part. Unbalanced MC's are seldom if ever mentioned in the literature, but if needed can be simulated by zero-padding in the appropriate places.

§4b. The C++ Standard Library contains a std::complex class. A candidate method for producing MC's might have been to form this sequence:

• ℂ₀ → double
• ℂ₁ → std::complex<double>
• ℂ₂ → std::complex<std::complex<double>>
• ℂ₃ → std::complex<std::complex<std::complex<double>>>
• et cetera

However, the C++ standard does not require that the std::complex class be compatible with any but three element types: float (single-precision floating point), double, and long double (extra-precision floating point). Thus std::complex is permitted to fail for anything else, even other numerical types. A typical implementation of std::complex is designed to be highly efficient and robust with only the three target numerical types, and it forfeits much applicability to other element types.

In particular, the version of std::complex available to the present author was that of LLVM, which is well-suited to its intended purpose, but which became useless when std::complex<std::complex<double>> was attempted for ℂ₂. In this regard, std::complex differs from genuine containers such as std::vector and std::list which nest without complaint.

§4c. Class multicomplex is nestable, because there was no reason to make it non-nestable; unclear is whether nestability useful in practice. Nested MC's of types of equal total depth, such as:

• using mcd5 = multicomplex<double,5>;
• using mcd32 = multicomplex<multicomplex<double,3>,2>;
• using mcd23 = multicomplex<multicomplex<double,2>,3>;

When nested MC's are used, view_k should be chosen for output, either directly or by way of operator<<. Functions view_i and view_j will not work, and even if modified would be confusing to use. This is because the subscript numbers of i and j would be unaware of the nesting, and each i_n or j_n symbol would be employed by different levels of nesting with different meanings.

The three types mcd5, mcd32, and mcd23 are likely to have the same layout in computer memory, and do so in the implementation available to the author (see layout_demo). Something else that will probably exhibit the same physical arrangement is a C-style array of the element type. Whether the C++ specification guarantees all these to have the same layout, however, is not clear. Constrast that class flat, which uses std::vector, will likely have a different layout because std::vector incurs overhead for changing the size of the data structure, even though flat never requests such a size change.

§4d. It is straightforward to adapt multicomplex to handle multi-quaternions. A quaternion has four elements rather than two, the multiplication rule is notably more complicated than that of complex numbers, and it must be remembered that multiplication is not commutative. Also, many objects in the code would need their names changed for clarity.

§5a. Complex numbers, hence multicomplex numbers, can be represented using ordinary matrices from linear algebra. Matrix addition or multiplication corresponds to MC addition or multiplication, and an MC has a reciprocal if the corresponding matrix has an inverse. Negating the matrix negates the MC, and transposing the matrix gives the conjugate of the MC.

The linear algebra model, although valid for computation, is not efficient in either execution time or memory usage. Still, it is sometimes helpful in answering theoretical questions, as linear algebra has developed a vast body of results.

In this report, the matrices are drawn by means of HTML tables, with some coloration to help show the structure.

§5b. In the ℂ₁ case, an MC named m can be written with either i or j notation:

equivalent:
m =	m0	+	m1i1
m =	m0j0	+	m1j1

and here is the matrix:

+ m0	− m1
+ m1	+ m0

Some authors prefer the following variant, which is equivalent, but it will not be pursued here:

+ m0	+ m1
− m1	+ m0

In the ℂ₂ case:

equivalent:
m =	m0	+	m1i1	+	m2i2	+	m3i2i1
m =	m0j0	+	m1j1	+	m2j2	+	m3j3

Of the ℂ₂ matrix, below:

• The upper left and lower right quadrants are the same as the ℂ₁ matrix.
• The lower left quadrant of the ℂ₂ matrix is like the upper left except every subscript has been increased by two.
• The upper right quadrant is the negative of the lower left.

+ m0	− m1	− m2	+ m3
+ m1	+ m0	− m3	− m2
+ m2	− m3	+ m0	− m1
+ m3	+ m2	+ m1	+ m0

The pattern continues with the ℂ₃ matrix, except that the subscripts in the lower left and upper right quadrants are now increased by four.

equivalent:

m =

m0

+

m1i1

+

m2i2

+

m3i2i1

+

m4i3

+

m5i3i1

+

m6i3i2

+

m7i3i2i1

m =

m0j0

+

m1j1

+

m2j2

+

m3j3

+

m4j4

+

m5j5

+

m6j6

+

m7j7

+ m0	− m1	− m2	+ m3	− m4	+ m5	+ m6	− m7
+ m1	+ m0	− m3	− m2	− m5	− m4	+ m7	+ m6
+ m2	− m3	+ m0	− m1	− m6	+ m7	− m4	+ m5
+ m3	+ m2	+ m1	+ m0	− m7	− m6	− m5	− m4
+ m4	− m5	− m6	+ m7	+ m0	− m1	− m2	+ m3
+ m5	+ m4	− m7	− m6	+ m1	+ m0	− m3	− m2
+ m6	− m7	+ m4	− m5	+ m2	− m3	+ m0	− m1
+ m7	+ m6	+ m5	+ m4	+ m3	+ m2	+ m1	+ m0

§5c. The scheme extends to ℂ_n, with the subscript increase each time being 2ⁿ⁻¹. There are formulas to calculate the number of positive and negative signs within the matrix:

	ℂ0	ℂ1	ℂ2	ℂ3	ℂ4	ℂ5	ℂn
number of positive signs	1	3	10	36	136	528	Pn = 22n−1 + 2n−1
number of negative signs	0	1	6	28	120	496	Nn = 22n−1 − 2n−1

As n becomes large, P_n ÷ N_n approaches unity.