Complex Numbers

[S-Erik]

The Idea of this Issue is to track the current progress on implementing complex numbers.

Following possibilites come to my mind:

• Introduce complex numbers as a type of this crate (struct Complex)
• Use the Complex type from the num_complex crate.

num_complex works fairly well with the ndarray crate, but for some simple operations such as computing the conjugate of an array you have to use the map function and apply the conj method to every element of the array. If we introduce complex numbers as a part of this crate we can implement functions like conj, dagger, etc. that operate on the whole array as well as on a single complex number, similar to the tranpose method t from ndarray.

What are your ideas on how to implement complex numbers?

@JErnestoMtz

[JErnestoMtz]

I think is worth creating our own Complex implementation, rapl does not have any dependencies and I feel that we should try to keep that way it if possible. Something like this should be fairly easy and effective:

struct Complex<T: Singed>{
    real: T,
    imag: T
}

impl <T: Singed> Complex<T> {
   fn conj(self)->Self{
    Complex { real: self.real, imag: -self.imag }
   }
   fn abs(self)->Self{
    todo!()
   }
   fn real(self)->Self{
    todo!()
   }
   fn imag(self)->Self{
    todo!()
   }
 }

Where the Signed trait would be something very similar to the num traits implementation, and all singed numeric types like i32, f32, etc. should implement that.

An other aspect is worth discussing about complex numbers, is if we should allow arithmetic between Complex<T> and non complex <T>, in the same way we allow arithmetic between Ndarr and scalars.

let a = Complex::new(2,3) + 1;
assert_eq!(a, Complex::new(3,3))

I think something like that is worth it. By default T: Signed would be real, but also we should add an extra trait to convert it to imaginary:

let a = Complex::new(2,3) + 1.imag();
assert_eq!(a, Complex::new(2,4))

[DeliciousHair]

Been thinking about this a bit, I think the signed bit may require some extra features. Since imag is not ordered over real, the signing part doesn't quite seem necessary from a conceptual point of view. More significantly though, if imag is signed (in the sense of u32 vs i32) then things as foundational as conjugation will become tricky at best, no?

[JErnestoMtz]

yes, I agree @DeliciousHair . In fact I find an implementation that I like a lot and solves that exact problem https://docs.rs/complex_algebra/0.1.8/complex_algebra/index.html

The basic complex type just required partial equality:

pub struct Complex<T: PartialEq>(pub T, pub T);

And all the other specific operations, like conj, abs, exp, trig functions and so on just require the minimal necessary traits; for example:

impl <T: Copy + PartialEq + Neg<Output = T>> Complex<T>{
   pub fn conj(&self)->C<T>{
        Complex(self.0, -self.1)
   }
}
impl<T: Copy + PartialEq + Add<Output = T> + Mul<Output = T>> Complex<T> {
    pub fn r_square(&self) -> T {
        self.0 * self.0 + self.1 * self.1
    }
}

I think doing it this way we can maintain generality, without sacrificing functionality. Although the create I mention goes a little bit to far, defining traits for cosine, sine, exponent, instead of just one Trig trait. So there are some small details that I want to change, but the basic idea is to keep the requirements in the functions not in the type

[DeliciousHair]

That looks good! Do you think it would also be worth implementing a way to house data in Euler notation? That is, a magnitude and phase angle directly in order to leverage various computations directly rather than the back and forth with trigonometry?

[S-Erik]

I tried to implement basics functions for our complex numbers, like new, conj, abs_squared, abs. The abs makes problems. My implementation of new, conj, abs_squared:

use std::ops::{Neg, Add, Mul};

#[derive(Debug, Clone, PartialEq)]
pub struct Complex<T: PartialEq>{
    pub real: T,
    pub imag: T
}

impl <T: PartialEq> Complex<T>{
    pub fn new(real: T, imag: T) -> Self {
        Self{real, imag}
    }
}

impl <T: Copy + PartialEq + Neg<Output = T>> Complex<T>{
    pub fn conj(&self)->Complex<T>{
         Complex::new(self.real, -self.imag)
    }
 }

 impl<T: Copy + PartialEq + Add<Output = T> + Mul<Output = T>> Complex<T> {

     pub fn abs_squared(&self) -> T {
         self.real * self.real + self.imag * self.imag
     }
 }

Now the abs function (does not compile):

 impl<T: PartialEq + Add<Output = T> + Mul<Output = T>> Complex<T> {
    pub fn abs(&self) -> T {
        (self.real * self.real + self.imag * self.imag).sqrt()
    }
 }

The issue is that the sqrt function is only implemented for the float types f32 and f64. Now one can implement the abs function for both types separately like this (this does compile):

impl Complex<f32> {
    pub fn abs(&self) -> f32 {
        (self.real * self.real + self.imag * self.imag).sqrt()
    }
}
impl Complex<f64> {
    pub fn abs(&self) -> f64{
        (self.real * self.real + self.imag * self.imag).sqrt()
    }
}

Personally, I do not like this solution because I suppose this will not be the only function where we encounter such a problem besides the obvious problem of code duplication.

But there is a solution. The num_traits crate defines a trait Float which basically means f32 or f64. So we can write something like this:

 impl<T: Float> Complex<T> {
    pub fn abs(&self) -> T {
        (self.real * self.real + self.imag * self.imag).sqrt()
    }
 }

Besides that, the num_trait crate defines other useful features such as nan or infinity. Similar the trait Pow of the num_traits crate is useful to calculate powers, which is a essential feature for rapl in my eyes.

Using num_traits involves having a dependency of course, but I honestly think that it is not a bad thing because we can build on established crates and existing features that are widely used without implementing these from scratch for rapl, although we of course depend on code from others. In a similar way I now think that we should use the num_complex crate for complex numbers and use it and modify parts to fit in rapl, e.g implementing the conj function on a complex Ndarr.

What do you guys think about this?

[JErnestoMtz]

I believe complex numbers is so commonly use that is worth making our own implementations that feels part of rapl in an almost native way; that includes having the same philosophy of simplicity and elegance for the API. There is always the possibility to use num_complex if someone does not like our implementation.

I created a temporal repository with my attempt of doing this: https://github.com/JErnestoMtz/rapl_complex

I haven't have time for writing all the tests, but in theory it already has most of the functionality num_complex, with a interface that I think fits better with rapl. I did use the Float trait of num_traits as you recommended @S-Erik.

Also, my math might be a bit wrong in some of the trig functions or exponents, I based it from some old notes form college.

I'm still debating about the complex struct itself, I use:

pub struct C<T: PartialEq>(pub T, pub T);

Because from a usability perspective, it is the fastest way to write a complex number C(1,2) . But I do see a point in using the originally proposed:

pub struct Complex<T: PartialEq>{
    pub real: T,
    pub imag: T
}

Even more because in rapl_complex you can do 1 + 2.im() which is also elegant, and would be independent of the struct.

What do you guys prefer?

[JErnestoMtz]

I created a new branch for the development of Complex numbers: https://github.com/JErnestoMtz/rapl/tree/complex

There are still a few things missing that I think must be finished before merging with main, the ones that I'm aware are:

• [x] complex sqrt.
• [x] Complex power (real to complex, complex to real, real to real, complex float to float) .
• [x] Display trait implementation.
• [ ] More rigorous testing.

Also noteworthy, I made some changes in the complex branch to the core rapl, mainly in the ops.rs file to allow arithmetic operations between Complex and Non-complex Ndarr without the need of explicitly casting types:

        let x = Ndarr::from([1, 2, 3]);
        let y = Ndarr::from([1.i(), 1.i(), 1.i()]);
        assert_eq!(&x + 1.i(), x + y);

[JErnestoMtz]

In commit 017a5eadf3ca657b8a4baa409ffca068080e0bf9 I just added sqrt I was worried about how to chose one of the two solutions of the polynomial, but apparently other implementations just chose the main branch every time.

I also created two helping functions inspired from the approx to facilitate testing with floating numbers.

[JErnestoMtz]

I'm closing this issue in favor of more specific issues.

I also decided to merge complex branch to main, because it had some considerable improvements over main in non complex related stuff. And the Complex implementation is at a point that some people might find it useful.