How can we create a Mathematica function that can accept 0 or more arguments?

How can we create a Mathematica function that can accept 0 or more arguments in the same way the Times function does?

Specially, I'd like to replace the Times function to work with Reals, Complexes, Quaternions or 4x4 matrices.

It looks like we can do the following:

MySpecialFunction[x1_, x2_, ..., xn_] := ... (* I don't know how to create this function yet *)

Unprotect[Times];
Times = MySpecialFunction;
Protect[Times];

True == (Times[1, 2, 3, 4, 5] == MySpecialFunction[1, 2, 3, 4, 5]);
True == (Times[2, 3] == MySpecialFunction[2, 3]);
True == (Times[] == MySpecialFunction[]);

If MySpecialFunction is, for example, Plus, it works.

Attempts

(* It looks like ## can help us to get all arguments, but I don't know how to use it to create the function. *)
In:= f[##] &[a, b, c, d]
Out= f[a, b, c, d]
• "Quaternions or 4x4 matrices" - Times[] is commutative, so you probably want to be careful there. May 6 '20 at 14:44
• Looks like a duplicate of this. May 6 '20 at 14:48
• @J.M., I agree. It is just that I'd like to unify the * (Times), ** (NonCommutativeMultiply) and . (Dot) operations at one. But I agree that I need to be careful. I'd like to apply the same concept to other functions as well. Like Sin[matrix] be actually MatrixFunction[Sin, matrix] instead, to simplify my analysis. May 6 '20 at 14:54
• For another approach, look up BlankSequence and BlankNullSequence for pattern matching the many arguments in the function's definition. May 6 '20 at 15:13

Intro

I assume in the following that we have a multiplication operation mytimes which is non-commutative, an addition operation Plus which is abelian and scalar multiplication Times.

We will want our operation to be distributive and linear this is the tricky thing in my eyes.

Let's code it

My favorite Mathematica book is Power Programming with Mathematica - The Kernel by David B. Wagner. You can find a nice example how one could code a linear operator there. This case is very similar.

If we have something like mytimes[c*somehead[z]] and if ourc is a scalar (any element of $$\mathbb{C}$$), then we want to move it out mytimes (homogeneity). We also want to preserve the order in case we have multiple arguments.

pattern1 = c1_*c2_ /; Element[c1, Complexes];
mytimes[y___] /; Length[Cases[{y}, pattern1]] != 0 := {y} /. {x1___, x2_*x3_, x4___} /; Element[x2, Complexes] :> x2*mytimes[x1, x3, x4];

The left hand side of the upper code ensures that we only apply the definition if we have at least one element where at least one scalar is multiplied with something else.

Next is distributivity:

pattern2 = c1_ + c2_;
mytimes[y___] /; Length[Cases[{y}, pattern2]] != 0 := {y} /. {x1___, x2_ + x3_, x4___} :> mytimes[x1, x2, x4] + mytimes[x1, x3, x4];

Analogous to the previous case we only want to match if at least one argument consists of a sum. If we have such an argument we let the pattern matching engine do the complicated stuff on the right hand side.

Lets see some stuff

With the upper definitions