I am programming commutative hypercomplex number class.
First of all: If that already exist, please tell. I have not found. PS: Quaternion exists but are not commutative.
I am trying to reproduce all that Matlab can do with real and complex, but with hypercomplex. The solution that I am working on is to define a class. If a more smart method may be used, please tell.
Many codes with complex numbers exist, the idea is to overload smartly to make working all existing program with Hypercomplex. This work might be used by many users.
please see at the end of the post the main maths principle found in Ablamowicz1996 _ Clifford Algebras with Numeric and Symbolic Computations
My Hypercomplex class works for basic operation: I have overloaded the mtimes(a,b) = a*b operation.
For instance, I can do:
A = MyClass(1);
B = MyClass(2);
A*B; % ->that works
As a main problem, for instance in 2D, I would like to be able to use
[A B C].*[D ; E]
while using the element-wise behavior of Matlab, and for any dimension , and with my definition of *. For instance, while dealing with matrices of compatible size but not identical:
A = [1 2 3]; B = [1;2]
Matlab gives
A.*B
ans =
1 2 3
2 4 6
That works also for N dimension when compatible.
How can I overload times(a,b) = a.*b , power, rdivide.. in order to have the element-wise behavior of Matlab?
Here is my Class
classdef Hypercomplex < handle
properties (SetAccess = public, GetAccess = public)
A
end
properties (SetAccess = private, GetAccess = private)
eps,nu,e1,e2
end
methods
function s = Hypercomplex(Z) %% Conctruc
% Z might be real, complex or hypercomplex
if isa(Z,'Hypercomplex') % is HyperComplex
s = Z;
elseif ~isreal(Z) && length(Z) == 1 % is Complex
s = Hypercomplex([real(Z) imag(Z) 0 0]);
elseif isreal(Z) && length(Z) < 5
s.A = zeros(1,4);
for ii = 1:length(Z) % is vecteur
s.A(ii) = Z(ii);
end
end
end
function s = epsnu(s)
s.eps = (s.A(1) - s.A(4)) 1i*(s.A(2) s.A(3));
s.nu = (s.A(1) s.A(4)) 1i*(s.A(2)-s.A(3));
end
function s = e1e2(s)
s.e1 = Hypercomplex([0.5 0 0 -0.5]);
s.e2 = Hypercomplex([0.5 0 0 0.5]);
end
%% Overload
function out = cos(a)
a.epsnu.e1e2;
out = cos(a.eps)*a.e1 cos(a.nu)*a.e2;
% I might also do for sin exp sqrt ...
end
function out = mtimes(a,b) % '*'
a = Hypercomplex(a).epsnu;
b = Hypercomplex(b).epsnu;
eps12 = 0.5*(a.eps * b.eps);% produit complex
nu12 = 0.5*(a.nu * b.nu );% produit complex
h1 = Hypercomplex([real(eps12) imag(eps12) imag(eps12) -real(eps12)]);
h2 = Hypercomplex([real(nu12) imag(nu12) -imag(nu12) real(nu12) ]);
out = h1 h2;
end
function out = mrdivide(a,b)
a = Hypercomplex(a);
b = Hypercomplex(b);
out = a * b.inverse;
end
function s = inverse(s)
s = Hypercomplex(s).epsnu.e1e2;
s = 1/(s.eps)*s.e1 1/(s.nu)*s.e2;
end
function out = plus(a,b)
a = Hypercomplex(a);
b = Hypercomplex(b);
out = Hypercomplex(a.A b.A);
end
function out = minus(a,b)
a = Hypercomplex(a);
b = Hypercomplex(b);
out = Hypercomplex(a.A - b.A);
end
end
end
Example i^2 = -1 :
>> Ni = Hypercomplex([0 1 0 0]) ;
>> Ni*Ni
Hypercomplex with properties:
A: [-1 0 0 0]
CodePudding user response：
I would think you need to loop over your arrays, and have a special case where, if for one dimension one of the inputs has size 1 and the other doesn't, then the same element for that array is used. A simple way to accomplish this is to index with (i-1)*step 1, where step is 0 if the array has size 1, and 1 otherwise.
One thing to note is that, in your overloaded method, indexing doesn't use your custom indexing definition (subsref and subsasgn), but the built-in one.
For example, assuming you only allow 2D arrays of your class:
function C = times(A, B)
szA = size(A);
szB = size(B);
if numel(szA)~=2 || numel(szB)~=2
error('Only 2D arrays are supported')
end
szC = max(szA, szB);
stepA = [matchSizes(szA(1), szC(1)), matchSizes(szA(2), szC(2))];
stepB = [matchSizes(szB(1), szC(1)), matchSizes(szB(2), szC(2))];
C = MyClass.empty(szC);
for jj = 1:szC(2)
for ii = 1:szC(1)
C(ii,jj) = A((ii-1)*stepA(1) 1,(jj-1)*stepA(2) 1) * B((ii-1)*stepB(1) 1,(jj-1)*stepB(2) 1);
end
end
% private helper function, keep in the same file
function step = matchSizes(szIn, szOut)
step = 1;
if szIn ~= szOut
if szIn == 1
step = 0;
else
error('Sizes not compatible')
end
end
NOTE: Code not tested.
NOTE2: This is one of the cases where 0-based indexing is more convenient!
CodePudding user response：
If MyClass has a property val to store numeric values I think (hope) that this should be possible (not tested):
function obj = times(A, B)
r = reshape([A.val], size(A)) .* reshape([B.val], size(B));
obj = arrayfun(@MyClass, r);
end
EDIT (because question is edited):
You can follow two ways to implement array of those numbers:
Array of objects
In this method each object just contains a hypercomplex number. It is the way that you have tried. As @CrisLuengo pointed out you can use loops to create the desired result or you can use a combination of vectorization and use of functions likearrayfunandcellfun. You need to know how to work with cell arrays and comma separated list and be familiar with verctorization operations like implicit expansion. Here I implementedplusandepsnuusing vectorization andcellfun. They can work with both scalar and object arrays:function out = plus(a,b) a = Hypercomplex(a); b = Hypercomplex(b); AA = reshape([a.A], [4, size(a)]); BB = reshape([b.A], [4, size(b)]]) CC = AA BB; A1 = reshape() out = cellfun(@Hypercomplex, num2cell(reshape(CC, 4, []).', 2)); end function s = epsnu(s) AA = reshape([s.A], 4, []); eps = num2cell(AA(1, :) - AA(4, :) 1i.*(AA(2, :) AA(3, :))); nu = num2cell(AA(1, :) AA(4, :) 1i.*(AA(2, :) - AA(3, :))); [s.eps] = eps{:}; [s.nu] = nu{:}; endObject containing array of numbers
Instead of definingAas a[4 x 1]array you can define each coordinate(x,y,z,et)as a separate property of class each can contain array of numbers. you need to overload thecatoperator to concatenate objects. As a result no object array is created but the underlying verctors are concatenated. You need tooverloadsubsrefandsubsasgnif you want to simulate element access. I think this method is simpler to implement than the first method.