[Libre-soc-dev] complex numbers as matrices

[Jacob Lifshay]

On Wed, Jun 30, 2021, 11:07 Richard Wilbur <richard.wilbur at gmail.com> wrote:

>
> 
> > On Jun 28, 2021, at 11:27, Luke Kenneth Casson Leighton <lkcl at lkcl.net>
> wrote:
> >
> > On 6/28/21, Richard Wilbur <richard.wilbur at gmail.com> wrote:
> >> o think about.  Can
> >> the vector engine handle multiplying a 3-dimensional matrix by a
> >> 4-dimensional matrix, yielding a 3-dimensional matrix result?
> >
> > yep.
> >
> > https://libre-soc.org/openpower/sv/remap/
> >
> > wait... *4* dimensional, no.
> >
> > only goes as high as 3 dimensions, XYZ.
> >
> > however... if you don't mind a manual for-loop on the 4th outer
> > dimension it's fine.
> >
> > i assume you didn't mean two 2-dimensional matrices, one 4x3 the other
> > 3x4 to create a 2D matrix 3x3 result.
> >
> > i believe you *actually* mean, "a matrix of NxNxNxN multiplied by one
> > of NxNxN" which makes no sense, but hey
>
> Let’s take a very simplistic example:  the circuit consists of a periodic
> voltage source, V(ω), connected across a purely inductive load, L.
> According to Ohm’s Law, V=IZ.  This is for the AC case a complex-valued
> phasor equation, V(ω)=I(ω)Z(ω), where Z(ω) = jωL is the impedance of the
> inductor.
>
> In this simplistic example, solving for the current, I(ω), in the
> inductive load is given by dividing both sides by the impedance, a complex
> number.
>
> If we treat complex numbers as scalars, this is a 1x1 matrix equation.
> (Hence, 0D since it is normal practice to ignore dimensions whose size is
> 1.)  If we represent complex numbers as 2x2 real matrices this becomes a
> 1x1x2x2 matrix equation, and thus 2D.
>
> When we consider a more complicated circuit with N elements and use
> Kirchoff’s Current or Voltage Laws to analyze it by nodes or loops, the
> result will be a set of N simultaneous equations in N variables.  The
> solution then involves an equation of
> [ 1xN ] = [ 1xN ] x [ NxN ]
> where every element is complex.  Hence for first-class complex treatment
> we have
> [ 1xNx1 ] = [ 1xNx1 ] x [ NxNx1 ]
> which yields under dimensional analysis
> 1D = 1D x 2D
> which is hardly a challenge for a vector machine unless N grows too large.
>
> On the other hand, if we choose to represent the complex numbers as 2x2
> matrices of reals, then we have the following
> [ 1xNx2x2] = [ 1xNx2x2 ] x [ NxNx2x2 ]
> which yields under dimensional analysis
> 3D = 3D x 4D
> and thus no longer necessarily so trivial, even for a vector machine.
>
> I would consider dedicating a mode to dealing with pairs of scalars that
> travel together as complex numbers.  The rules are slightly different than
> regular [ 1x2 ] vectors but they are very important in science and
> engineering and, this feeds into a different conversation on this mailing
> list, the normal result of a Fourier transform.
>
> With a mode to treat register pairs as the real and imaginary parts of
> complex numbers we can have double-precision complex number support for
> free.  Multiplication:
> complex x complex = 4 scalar products, 1 negation, 2 additions
> complex x real (or complex with imaginary part =0) = 2 scalar products
> (this is already covered by your vector multiply with one factor a scalar)
> complex x imaginary (or complex with real part =0) = 2 scalar products, 1
> negation, destination swap
> complex + complex = 2 scalar additions (covered by vector addition)
> complex - complex = 2 scalar negations, 2 scalar additions (covered by
> vector operations)
>

wouldn't this just be a vec2 subtraction? no need for separate negation and
addition...

abs(complex= a + jb) = sqrt(complex x conjugate(complex))
>   = sqrt((a + jb)(a - jb)) = sqrt(a^2 + b^2)
> (Same as vector magnitude.)
>
> argument(complex= x + jy) = arctan(y/x) (this can have a meaningful value
> for all real values of x,y except (0,0), same as 2D vector)(I recommend
> consulting the implementation of atan2)[*]
> Normally this is defined with a principal branch covering the range θ in
> (-π, π].
>
> conjugate(complex) = complex* = (a + jb)* = a - jb (entails negation of
> imaginary part)
>
> reciprocal(complex) = (1/abs(complex)) x conjugate(complex)
>

that's wrong, the correct expression is:
z = re+i*im
reciprocal(z) = conjugate(z) / (re^2 + im^2)

>
> equality (for z1 = a1 + jb1, z2 = a2 + jb2)
> z1 == z2 <=> ((a1 == a2) and (b1 == b2)) (this is the same as vectors)
>
> ordering is more complicated than normal scalars:  Two fairly simple
> possible schemas come to mind and both involve the polar form.  z1 = (r1,
> θ1), z2 = (r2, θ2)
> r-major:  z1 < z2 <=> ((r1 < r2) or ((r1 == r2) and (θ1 < θ2))),
> θ-major:  z1 < z2 <=> ((θ1 < θ2) or ((θ1 == θ2) and (r1 < r2)))
>

I would argue that, if we provide ordering for complex numbers, we should
use real-major form since it's much simpler to compute and since complex
numbers don't normally have an ordering defined.
z1 = r1 + i*i1, z2 = r2 + i*i2
z1 < z2 <=> (r1 < r2) || (r1 == r2 && i1 < i2)

This is the same order you would get in Rust for:
#[derive(PartialEq, PartialOrd)]
pub struct MyComplex {
    real: f32,
    imag: f32,
}
or:
pub type MyComplex = (f32, f32);

or in C++ for:
struct MyComplex {
    float real, imag;
    auto operator <=>(const MyComplex &) const = default;
};

Jacob