Defining a non-standard algebraic number

Would this work? I am using a special symbol for your $c$, to help distinguish it from any other variable:

Clear[\[ConstantC]] \[ConstantC] /: \[ConstantC] \[ConstantC] = 1 \[ConstantC] /: \[ConstantC] I = -I \[ConstantC] 

This looks like the following in a notebook:

With these definitions, I can evaluate:

and obtain:

You might try embedding into a larger matrix in a way that the relations work out. So the imaginary unit becomes {{0,1},{1,0}}, real values become diagonal 2x2 matrices, and for this particular c it looks like {{1,0},{0,-1}} will work. This will involve simple replacements, some array flattening and unflattening, and a simple undoing of the replacements. Also important is to realize that we really have four generators, not 3, insofar as resulting 2x2 matrices will need to be reconstructed as a real part, an imaginary part, a c part, and one other reverse diagonal which can be in terms of either "I*c" or "c*I". Since multiplications is commutative we cannot use these explicitly and so I will just use a surrogate variable, called ci, for the latter. Also I will use i from the beginning to represent out imaginary unit, rather than Complex[0,1]. So here is how we proceed.

Start with replacements that take us from scalar to 2x2 matrix form, and reverse.

replacements = {i -> {{0, 1}, {1, 0}}, c -> {{1, 0}, {0, -1}}}; cirForm[mat_?MatrixQ /; Dimensions[mat] === {2, 2}] := (mat[[1, 1]] + mat[[2, 2]])/2 + c*(mat[[1, 1]] - mat[[2, 2]])/2 + i*(mat[[1, 2]] + mat[[2, 1]])/2 + ci*(mat[[1, 2]] - mat[[2, 1]])/2 

We create an example using two 4x4 matrices containing real, imaginary, and c parts.

SeedRandom[1234]; rand1 = RandomInteger[{-10, 10}, {3, 4, 4}]; rand2 = RandomInteger[{-10, 10}, {3, 4, 4}]; mat1 = {1, i, c}.rand1 mat2 = {1, i, c}.rand2 (* Out[1606]= {{-9 - 6 c + 6 i, 10 - 8 c + 8 i, 10 + 6 c + 10 i, -9 - 6 c - 3 i}, {-10 + 6 c + 4 i, 6 - 9 c + 6 i, 7 + 3 c - 4 i, -6 + 7 c - 9 i}, {-5 - c + 5 i, 9 + 10 c - 3 i, 5 + c - 2 i, 6 - 4 c + 3 i}, {9 + 4 c + 3 i, -8 + 8 c - 9 i, 2 + c + 5 i, -5 - 8 c - 2 i}} Out[1607]= {{-6 - 8 c - 5 i, 6 - 6 c - i, 5 - 8 c + 6 i, 1 + 10 i}, {-1 - 2 i, 4 + 6 c - 10 i, 2 + 4 c + 7 i, 9 - 6 c - i}, {9 + 2 c + 8 i, 10 + c - 9 i, 1 + 9 c - 8 i, -9 - 9 c + 4 i}, {-2 + 4 c - 5 i, 7 + 8 i, -8 + 4 c - 2 i, -8 + c + 2 i}} *) 

Now flatten into 8x8s.

mat1Full = ArrayFlatten[mat1 /. replacements, 2] mat2Full = ArrayFlatten[mat2 /. replacements, 2] (* Out[1608]= {{-15, -3, 2, 18, 16, 20, -15, -12}, {-3, -3, 18, 18, 20, 4, -12, -3}, {-4, -6, -3, 12, 10, 3, 1, -15}, {-6, -16, 12, 15, 3, 4, -15, -13}, {-6, 0, 19, 6, 6, 3, 2, 9}, {0, -4, 6, -1, 3, 4, 9, 10}, {13, 12, 0, -17, 3, 7, -13, -7}, {12, 5, -17, -16, 7, 1, -7, 3}} Out[1609]= {{-14, -11, 0, 5, -3, 11, 1, 11}, {-11, 2, 5, 12, 11, 13, 11, 1}, {-1, -3, 10, -6, 6, 9, 3, 8}, {-3, -1, -6, -2, 9, -2, 8, 15}, {11, 17, 11, 1, 10, -7, -18, -5}, {17, 7, 1, 9, -7, -8, -5, 0}, {2, -7, 7, 15, -4, -10, -7, -6}, {-7, -6, 15, 7, -10, -12, -6, -9}} *) 

Take the product and repartition into a 4x4 comprised of 2x2s.

matProductFull = mat1Full.mat2Full (* Out[1600]= {{757, 724, -192, -272, 386, -200, -109, 236}, {288, 425, 152, -340, 496, 38, -116, 377}, {357, 303, -237, -151, 261, -97, -95, 185}, {365, 245, -313, -601, 241, 57, 83, 431}, {105, 58, 372, -30, 127, -101, -92, 53}, {90, -69, 296, 162, -151, -259, -231, -130}, {-88, 131, 6, 65, 43, 470, 53, 26}, {-99, 102, 25, 186, -166, 53, -212, -259}} *) matProduct = Partition[matProductFull, {2, 2}] (* Out[1601]= {{{{757, 724}, {288, 425}}, {{-192, -272}, {152, -340}}, {{386, -200}, {496, 38}}, {{-109, 236}, {-116, 377}}}, {{{357, 303}, {365, 245}}, {{-237, -151}, {-313, -601}}, {{261, -97}, {241, 57}}, {{-95, 185}, {83, 431}}}, {{{105, 58}, {90, -69}}, {{372, -30}, {296, 162}}, {{127, -101}, {-151, -259}}, {{-92, 53}, {-231, -130}}}, {{{-88, 131}, {-99, 102}}, {{6, 65}, {25, 186}}, {{43, 470}, {-166, 53}}, {{53, 26}, {-212, -259}}}} *) 

Rewrite in "cir" form.

cirmatProduct = Map[cirForm, matProduct, {2}] (* Out[1602]= {{591 + 166 c + 218 ci + 506 i, -266 + 74 c - 212 ci - 60 i, 212 + 174 c - 348 ci + 148 i, 134 - 243 c + 176 ci + 60 i}, {301 + 56 c - 31 ci + 334 i, -419 + 182 c + 81 ci - 232 i, 159 + 102 c - 169 ci + 72 i, 168 - 263 c + 51 ci + 134 i}, {18 + 87 c - 16 ci + 74 i, 267 + 105 c - 163 ci + 133 i, -66 + 193 c + 25 ci - 126 i, -111 + 19 c + 142 ci - 89 i}, {7 - 95 c + 115 ci + 16 i, 96 - 90 c + 20 ci + 45 i, 48 - 5 c + 318 ci + 152 i, -103 + 156 c + 119 ci - 93 i}} *) 

Your system is called "split-quaternions", or simply an algebra of $2\times2$ matrices.

You can use the following code to realize this system. In this code I denoted your constant c as j.

Clear["Global`*"]; Unprotect[Power]; Power[0, 0] = 1; Protect[Power]; $PrePrint =.; FormatValues@Matrix = {}; Hold[MakeBoxes[Matrix[a_], StandardForm] ^:= Block[{Internal`$ConvertForms = {}}, TemplateBox[{MakeBoxes[MatrixForm@a][[1, 1, 3]]}, "Matrix", DisplayFunction -> (RowBox@{style@"(", "\[NoBreak]", #, "\[NoBreak]", style@")"} &)]]] /. style[string_] -> StyleBox[string, FontColor -> Red] // ReleaseHold Matrix /: Plus[Matrix[A_?SquareMatrixQ], Matrix[B_?SquareMatrixQ]] := Matrix[KroneckerProduct[A, IdentityMatrix[LCM[Length[A], Length[B]]/Length[A]]] + KroneckerProduct[B, IdentityMatrix[LCM[Length[A], Length[B]]/Length[B]]]] Matrix /: Dot[Matrix[A_?SquareMatrixQ], Matrix[B_?SquareMatrixQ]] := Matrix[KroneckerProduct[A, IdentityMatrix[LCM[Length[A], Length[B]]/Length[A]]] . KroneckerProduct[B, IdentityMatrix[LCM[Length[A], Length[B]]/Length[B]]]] Matrix /: Matrix[x_?MatrixQ] := First[First[x]] /; x == First[First[x]] IdentityMatrix[Length[x]]; Matrix /: x_?NumericQ + Matrix[y_] := Matrix[x IdentityMatrix[Length[y]] + y]; Matrix /: x_?NumericQ Matrix[y_] := Matrix[x y]; Matrix /: Matrix[x_]*Matrix[y_] := Matrix[x . y] /; x . y == y . x; Matrix /: Power[Matrix[x_?MatrixQ], 0] := 1; Matrix /: Re[Matrix[x_]] := Matrix[Re[x]]; Matrix /: Power[Matrix[x_?MatrixQ], y_?NumericQ] := Matrix[MatrixPower[x, y]]; Matrix /: Power[Matrix[x_?MatrixQ], Matrix[y_?MatrixQ]] := Exp[Log[Matrix[x]] . Matrix[y]]; Matrix /: Log[Matrix[x_?MatrixQ], Matrix[y_?MatrixQ]] := Log[Matrix[x]]^-1 . Log[Matrix[y]] Matrix /: Scalar[Matrix[x_?MatrixQ]] := Tr[x]/Length[x] Scalar[x_?NumericQ] := Re[x] $Post2 = FullSimplify[# /. j -> Matrix[( { {-1, 0}, {0, 1} } )] /. i -> Matrix[( { {0, -1}, {1, 0} } ) ] /. f_[args1___?NumericQ, Matrix[mat_], args2___?NumericQ] :> Matrix[MatrixFunction[f[args1, #, args2] &, mat]] /. Matrix[{{a_, c_}, {d_, b_}}] :> (a + b)/ 2 + (b - a)/2 j + (c - d)/2 i - (c + d)/2 i . j] &; $Post = Nest[$Post2, #, 3] &; 

Since it is not commutative, you have to use dot product . to multiply the elements. All standard functions are applicable.

In:=Log[i + j + 1] Out: -i + j In:=Exp[i + j] Out:=1 - i + j In:=i^j Out:=Cosh[\[Pi]/2] + i . j Sinh[\[Pi]/2] 

For non-commutative multiply, we can use NonCommutativeMultiply, or ** to define the algebra:

Clear[c]; c /: c ** c := 1 c /: c ** (a_ + b_) := c ** a + c ** b c /: c ** (a_Real b_) := a c ** b c /: c ** a_Complex := Conjugate[a] ** c rules = {a_Complex ** c + b_Complex ** c :> (a + b) ** c, c + b_Complex ** c :> (1 + b) ** c}; (* any other properties of c *) (c + I) ** (2 + 3 I + c) ** c //. rules (*(2 - 2 I) + (-2 + 2 I) ** c*)