Need common trait for floating and complex numbers

[cuviper]

From @carrutstick on July 18, 2017 23:14

This issue in in the vein of issue #209, but a bit more ambitious; hopefully I am not too far out of scope for this project.

Any mathematical function that accepts a real number can technically also accept a complex number, but currently there is no trait we can use to define a function that is generic over real and complex numbers (unless I'm missing something).

I see two possible solutions:

1. Shoehorn complex numbers into the existing Float trait (minimal work)
2. Refactor the float trait so that we separate the parts that refer to mathematical objects (the trig functions, for instance) from the parts that assume specific types of implementation (min_value(), max(), etc.).

Option 1

Certain functions in the Float trait may be ambiguous for complex numbers, such as the min_value() function. Broadly, there are three classes of functions in the trait:

1. |Self| -> Self
2. || -> Self
3. |Self| -> bool

Obviously the trigonometric functions of class 1 are not ambiguous (once we pick a branch cut), but functions like floor() and frac() are less obvious. I think a good rule of thumb for such functions is to have them apply independently to each part of the complex number, so that

Complex::new(x, y).floor() == Complex::new(x.floor(), y.floor())

For class 2, I'm much less certain, but I think that these should generally produce values that are mathematically equivalent to their real counterparts, so that

Complex<F>::min_value() == Complex<F>::new(F::MIN, F::zero())
Complex<F>::infinity() == Complex<F>::new(F::infinity(), F::zero())

For class 3 I don't think there's a good rule of thumb, and we would probably need to invent rules on a case-by-case basis. For instance, we would probably want a complex number to have is_nan() == true if either part was nan, but what if the other part is infinite? Do we also have is_infinite() == true? Or do we decide that something that is nan cannot be infinite?

There is also the special case of the min and max functions. I think the only consistent way to handle these is to just have them act on the real line, like class 2.

Option 2

Maybe someone who knows more abstract algebra than me can come up with better names, but we could have traits similar to Haskell's system, where RealFloat describes any kind of binary floating point number, whereas Floating is a more abstract description of any object suitable for use with trigonometric functions, etc. (this would include complex numbers, auto-differentiation primitives, diagonalizable matrices, etc.)

Obviously this option would be more work, and would result in at least a few breaking changes, but it would also be great to have the ability to program more generically with mathematical objects.

Copied from original issue: rust-num/num#321

[cuviper]

From @termoshtt on July 19, 2017 5:13

I want the common trait in num-traits crate, too.

FYI, I take another option in my crate ndarray-linalg for somewhat different objective: I define a trait Scalar to capture both real and complex numbers to define exp() for both types.

fn some_math<A: Scalar>(a: A) -> A {
  a * a.exp()
}

This function can be called with f64, f632, Complex<f64>, and Complex<f32>. Large part of math functions can be defined naturally for both real and complex numbers, e.g. sin, exp, ln, abs, and so on, and I think this should be involved in num-traits.

[cuviper]

There was a similar discussion on the forums. I don't recall anything coming out of this though. https://users.rust-lang.org/t/floating-point-traits-roundtable/5382

I would lean more towards your Option 1, if only because it's more conservative¹. We've already made some of the choices you're pondering in impl<T: Clone + Float> Complex<T>, so I would just forward Float for Complex<T> to those. I'd almost say to just move those current methods to the Float implementation, but inherent methods get special precedence so I think that would be a breaking change. Your idea to favor reals for "class 2" sounds fine to me. I really don't know what to do with integer_decode though...

I see the value in the more ambitious Option 2, but I think that would be better explored in a new crate of your own, if you're so inclined.

¹ (If crates could have traits, num would definitely impl Conservative.)

[cuviper]

From @carrutstick on August 10, 2017 14:40

Hi @cuviper, thanks for the input. I've made it as far as a Num implementation for Complex, but now I'm stuck on the PartialOrd requirement for Float. I'm not sure that there is a way to define partial_cmp for Complex (or at least not one that agrees with PartialEq).

Any thoughts on this? Is the PartialOrd constraint actually used anywhere?

[cuviper]

From @Libbum on August 10, 2017 14:45

Mathematically you should take the norm(), which would give you Complex<T> -> T, then you can sort that value (since any possible T will already implement PartialOrd). I think ... someone else may need to confirm that.

[cuviper]

From @carrutstick on August 10, 2017 14:50

@Libbum, the problem with that is that partial_cmp needs to return one of Less, Greater, or Equal, and the Equal result is inappropriate for two complex numbers that just happen to have the same norm. Maybe it would be better to just have the comparison return None in a large set of cases (such as when the norm is equal, but the numbers are not), but I am worried that this would lead to a lot of unexpected errors.

[cuviper]

From @Libbum on August 10, 2017 14:57

Ahh yes, valid point. I'll think about it a bit more; but hopefully the others have some suggestions too. I'd definitely agree that returning a large set of possible Nones would be worrisome though.

[cuviper]

Using norm() means that -1 > 0, but I think we'd at least want to be consistent with real comparisons. A linear ordering is a little better - for a+bi and c+di, first compare a and c, and if equal compare b and d. But I think you can still derive contradictions from this.

Note that C++ and Python at least do not provide ordering for complex numbers. It's a breaking change to remove PartialOrd from Float though.

[cuviper]

From @Libbum on August 11, 2017 9:20

OK, so I've done a bit more digging here.

I think this discussion on the mathematics stack exchange summarises the initial questions we have fairly well.

1. Complex numbers have no natural (linear) ordering that is compatible with an ordered field structure.
2. One can impose total order on a complex list via a lexicographical order.

But what are the pitfalls of a lexical ordering? The contradictions that @cuviper speaks of come into play when one wants to multiply complex values.

For real values: If a < b and 0 < c then ac < bc (1) Note that i may be written as 0 + i1 while 0 = 0 + i0. Therefore, by definition, 0 < i. Assume (1) holds for complex numbers and the lexicographic order. Obviously -1 < 1. Multiplying by i we get -i < i. Multiplying by i one more time produces -i² < i². By definition, this implies 1 < (-1) which is absurd.

(Further discussion here)

Most people seem to have the feeling that although it's possible to make an order on the complex numbers, such an order would not be very useful. Whether not very useful can turn to dangerous in the hands of someone not aware of this should at least be considered.

---

So the question we need to answer is are we OK with implementing this? I've taken a look at sorting in Matlab and other tools and they all seem pretty sketchy when it comes to their choices of ordering - which is now obvious why. But they do end up implementing some form of lexical ordering for complex values even if they contradict themselves:

Matlab> sort ([-4+i -3+i])
-3+i -4+i 
Matlab> -4+i < -3+i
true

From where I sit there are two choices:

1. Implement lexical ordering as @cuviper suggests, but stress in the documentation that this ordering has pitfalls.
2. Don't implement PartialOrd for the complex type and start investigating an alternate route to treat real and complex values together in a generic manner.

[cuviper]

Even though python itself doesn't order complex numbers, numpy uses lexical order:

numpy.sort_complex:

Sort a complex array using the real part first, then the imaginary part.

numpy.sort:

The sort order for complex numbers is lexicographic. If both the real and imaginary parts are non-nan then the order is determined by the real parts except when they are equal, in which case the order is determined by the imaginary parts.

And numpy also uses this order in element-wise comparisons of complex arrays.

I think that's enough precedent for me, but let's indeed mention the pitfalls in documentation.

[aplund]

I think the approach taken in nalgebra is the most sane one. Shoehorning complex numbers into single floats cannot sanely capture the full structure of complex numbers. For example, how are you meant to represent complex conjugation with a float type?

It seems that in the discussion here there has been some confusion about what is and isn't common between complex and real numbers. Complex number cannot be totally ordered (as has been demonstrated above). So total ordering should be a separate trait that applies to reals but not complex numbers.

The issue about functions is a little bit subtle too. What I presume people are talking about is standard field operations (addition, subtraction, multiplication and division). These are shared between complex and reals. But functions in general have completely different properties between the two types of number (see https://en.wikipedia.org/wiki/Complex_logarithm). So there should also be a division along the lines of a "Field" trait and some kind of "Single valued functions" trait which can be used for both.

The sanest approach would be to take time to filter out what specific traits are desired are and separate out the functionality accordingly.

Also, I think Complex itself should be a trait rather than a concrete type given that Floats (i.e. Reals) and Integers are traits.

[ReidAtcheson]

I think the approach taken in nalgebra is the most sane one. Shoehorning complex numbers into single floats cannot sanely capture the full structure of complex numbers. For example, how are you meant to represent complex conjugation with a float type?

Complex conjugate of a real number is a real number. It would be a no-op.

It seems that in the discussion here there has been some confusion about what is and isn't common between complex and real numbers. Complex number cannot be totally ordered (as has been demonstrated above). So total ordering should be a separate trait that applies to reals but not complex numbers.

Agree. For people (like me) wanting a common trait for float & complex it's to write better generic functions when it makes sense to do so. That way I don't have to copy+paste a function into two copies one for real one for complex even though they are in all other respects identical. Naturally, if the ordering properties of real numbers is critical to an algorithm then no trait will magically make it work for complex numbers. There are many algorithms that do not rely on ordering of reals to work, though.

The issue about functions is a little bit subtle too. What I presume people are talking about is standard field operations (addition, subtraction, multiplication and division). These are shared between complex and reals. But functions in general have completely different properties between the two types of number (see https://en.wikipedia.org/wiki/Complex_logarithm). So there should also be a division along the lines of a "Field" trait and some kind of "Single valued functions" trait which can be used for both.

Could you be more specific? Almost every special function defined on both complex & real number has the property that the real version of the special function is equivalent to the complex version evaluated on a complex number with imaginary part zero. You have to specify branch cuts for some special functions to be well-defined on complex field, but that really isn't that problematic.

Real & Complex numbers both have very similar algebraic properties, and the completion of those algebraic properties to transcendental functions results in the same special functions (sin,cos,sqrt,log,...) up to branch cuts. using only these algebraic properties one can express many useful algorithms - in my case linear algebra but that is not the only domain where this is true.

Aside from mathematical properties complex numbers are usually implemented with floating point components, therefore it's also very useful to be able to query floating point properties (e.g. like how we can do f32::epsilon()) for things like error control.