au icon indicating copy to clipboard operation
au copied to clipboard

Add Multiplication Support for Complex Types

Open MatthiasJeuthe opened this issue 11 months ago • 9 comments

Thank you very much for sharing your code.

The multiplication operator in Quantity is only enabled if T suffices std::is_arithmetic. However, my use case involves the use of complex numbers, e.g. std::complex. What about allowing such types as well?

As a quick work-around I adapted the code as follows:

    template<typename T, typename = void>
struct has_real_and_imag : std::false_type
{
};
template<typename T>
struct has_real_and_imag<T, std::void_t<decltype(std::declval<T>().real()), decltype(std::declval<T>().imag())>> : std::true_type
{
};
template<typename T>
struct is_complex : has_real_and_imag<T>
{
};

// Scalar and Complex multiplication.
template<typename T, typename = std::enable_if_t < std::is_arithmetic<T>::value || is_complex<T>::value >>
friend constexpr auto operator*(Quantity a, T s)
{
	return make_quantity<UnitT>(a.value_ * s);
}

MatthiasJeuthe avatar Mar 18 '24 20:03 MatthiasJeuthe

Thanks for bringing this to our attention! Sorry I didn't get to it earlier, but this week has been really busy at work and at home. I will play around with this as soon as I get a chance --- possibly this weekend. I'll update the ticket once I understand the prospects for a solution.

chiphogg avatar Mar 20 '24 12:03 chiphogg

Thank you very much for taking care of this issue. From what I can see now, adaptions must involve:

  • Adding support for all arithmetic operations with complex-valued types.
  • Adding real() and imag() functions to complex-valued quantities. They should return the T::value_type with the same unit as the Quantity.

Maybe it helps.

MatthiasJeuthe avatar Mar 20 '24 19:03 MatthiasJeuthe

I looked into this some over the weekend. It turns out to be surprisingly tricky and subtle.

First, why are these std::is_arithmetic SFINAE instances here in the first place? It's to prevent ambiguous overloads of operator* and friends when we multiply two Quantity types. The intent wasn't so much to restrict scalars to only arithmetic types. Rather, the goal was more like "any valid rep", but unfortunately, we don't have a good trait or concept to express that constraint (see #52).

Although I don't know how to express what we do want, there are a bunch of types we definitely don't want:

  • Any Quantity or QuantityPoint
  • Any unit
  • Magnitudes
  • etc.

I thought we might be able to get away with switching to a "deny list" approach that excludes these types, using "any empty type" as a proxy for our monovalue types. Unfortunately, this doesn't work: if users define a custom type with multiplication or division operators which work with Quantity types, we'll get an ambiguous overload, because the Quantity version of the operator will think that custom type is OK to use as a rep.

Next time I take a crack at this, I'll probably try a simpler approach that keeps the "allow list", but simply expands it to std::complex on arithmetic types. I'm guessing that would unblock your main use case, right? My next several weekends are pretty packed, so I think it'd be about a month before I can pick up this thread.

As for adding real() and imag() member functions, I don't think this would be appropriate. It would add too much complexity to the Quantity interface for a somewhat niche use case. I'd recommend a free function approach, where you write an API something like real(x) that can be made to work both for any x that has .real(), and for Quantity<U, std::complex<T>>.

chiphogg avatar Mar 24 '24 21:03 chiphogg

Okay, thank you very much for trying to solve it and your explanations. Adding std::complex will not resolve my issue, since I am also using complex types from other libraries. However, they all share the same real() and imag() functions, that's why I proposed a template specialization covering this.

I am not sure whether it is worth to add it to the Quantity inferface. At least for me, this is a quite ordinary use-case😄

MatthiasJeuthe avatar Mar 25 '24 12:03 MatthiasJeuthe

That's really helpful, and I'm grateful that you posted it before I went and implemented something that wouldn't have solved your use case!

I'll mull this over and see what strategies my subconscious can come up with. Maybe the fact that I'm so swamped for the next several weeks will be a blessing in disguise, giving my subconscious more opportunity to digest the problem. :slightly_smiling_face:

chiphogg avatar Mar 27 '24 01:03 chiphogg

Update: everyone here got sick, so we had to cancel this weekend's plans. But the silver lining is that it gave me the chance to try out my ideas!

Try the code in #229. It passes all of the new and existing unit tests, and all of the tests internally in Aurora's repo. Additionally, I verified that it has no measurable impact on compile time.

If it fixes your problem, then I'll take it out of draft mode and get it reviewed!

chiphogg avatar Mar 30 '24 21:03 chiphogg

It seems to work now for std::complex, very good. However, I face another issue (unrelated to this one) when I try to use thrust::complex from the CUDA library.

au::Quantity<au::Thickness, Type> a = au::meters(thrust::complex(1.0, 2.0));

Using MSVC and the commit you pointed out, I receive the following compilation error: "au\magnitude.hh(525,34): error C3615: Die constexpr-Funktion "au::detail::get_value_result" kann keinen konstanten Ausdruck ergeben"

MatthiasJeuthe avatar Apr 01 '24 07:04 MatthiasJeuthe

Interesting!

What is au::Thickness? I can't find it in the library, so I assume it's a custom unit you added. (If so, I'd recommend giving it a unit-name rather than a dimension-name like Thickness. The library is designed to read as though whatever goes in that slot is a unit.)

Anyway, the error itself was very suggestive. get_value_result() has to do with unit conversions. It turns out there's a whole lot here that just doesn't work well yet with anything that isn't std::is_arithmetic types. For example:

constexpr auto length = meters(std::complex<double>{1.0, 2.0});
const auto length_m = length.in(meters);          // Passes
const auto length_mm = length.in(milli(meters));  // Fails horribly

I'm going to need to figure out how to approach this problem in general... :thinking:

chiphogg avatar Apr 21 '24 02:04 chiphogg

@MatthiasJeuthe, could you please check the draft PR #278, and see if it fixes your problem? If it does, I'll get it ready for more thorough review.

chiphogg avatar Aug 04 '24 12:08 chiphogg

Thank you very much for working on the issue. I am quite busy at the moment. I try to test it during the weekend.

MatthiasJeuthe avatar Aug 13 '24 20:08 MatthiasJeuthe

Multiplication and division by a scalar appear to be functioning correctly; I verified this using both std::complex and cuda::std::complex. However, multiplying by a complex-valued quantity isn't supported yet. I assume that wasn't the goal of this fix, correct?

For example:

au::Quantity<au::Meters, cuda::std::complex<double>> test = au::meters(cuda::std::complex<double>{1.0, 2.0}); test *= au::meters(cuda::std::complex<double>{2.0, 2.0});

This results in the following error: C2338 static_assert failed: 'We don't support compound multiplication of integral types by floating point'.

MatthiasJeuthe avatar Aug 17 '24 06:08 MatthiasJeuthe

I hadn't considered that test case, but I think we should include it in this PR.

One small nitpick. As written, the test case could never work, because you can't pass au::meters(...) as the argument to *= for any type. After all, the left hand side must be the same type before and after, but multiplying by meters will always change the unit, and thus also the type.

That said, I've added some tests that multiply by a "raw" complex number. Your feedback was really helpful in uncovering an issue with complex support, and I think I've fixed it now. Try the latest on #278!

chiphogg avatar Aug 17 '24 15:08 chiphogg

@MatthiasJeuthe, let me know if you can come up with any more test cases we should include. Otherwise, I'll fix this up for review, likely sometime soon after my CppCon talk next Wednesday. :slightly_smiling_face:

chiphogg avatar Sep 10 '24 19:09 chiphogg

@chiphogg I believe you're ready to proceed with your review. Best of luck with your talk—I’ll definitely watch it on YouTube!

MatthiasJeuthe avatar Sep 10 '24 20:09 MatthiasJeuthe