We need a universal method of constructing Complex type from Real.

[cosurgi]

@ckormanyos mentioned that math does not depend on multiprecision. So we cannot simply #include <boost/multiprecision/complex_adaptor.hpp> and go ahead. Also there is a multitude of real types available within Boost. See for example this short program:

#include <boost/config.hpp>
#include <boost/cstdfloat.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/mpfi.hpp>
#include <boost/core/demangle.hpp>

template<typename T>
void print_sin_1() {
    using std::sin;
    std::cout << sin(T(1)) << " type: " << boost::core::demangle(typeid(T).name()) << std::endl;
}

template <typename... T>
constexpr inline auto for_each = [](auto&& f) {
    (f.template operator()<T>(), ...);
};

int main() {
    auto for_each_test_type = for_each<
              float
            , double
            , long double
            , boost::multiprecision::float128
            , boost::multiprecision::cpp_bin_float_50
            , boost::multiprecision::cpp_bin_float_quad
            , boost::multiprecision::cpp_dec_float_100
            , boost::multiprecision::number<boost::multiprecision::cpp_bin_float<200> >
            , boost::multiprecision::mpfr_float_100
            , boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<200> >
            , boost::multiprecision::mpfi_float_100
            , boost::multiprecision::number<boost::multiprecision::mpfi_float_backend<200> >
        >;

    for_each_test_type([]<typename T>() {
        print_sin_1<T>();
    });
}

Also this code snippet should be nice to test the solution to the problem which we seek here.

Not mentioning the types in the works such as double_float and quad_float (I believe complex_adaptor is for them).

Now we need something which, in a universal way, constructs a Complex type for each of them. We have a hackish tool for that here already used in the code in several places.

@Lagrang3 already noted in the comments that it is needed.

[cosurgi]

@ckormanyos I suspect that the details of the solution should be defined inside Boost.Multiprecision folder. In a way I see it as Boost.Multiprecision does not provide enough type traits for the Complex types available in there.

The user would include for example boost/multiprecision/mpc.hpp before FFT, to get a specialization of a class which solves this issue. On FFT side in Boost.Math we will only use an empty template without any specialization. User has to include boost/multiprecision/mpc.hpp (or any other Complex type) before calling FFT.

Does it mean that we need to fork multiprecision here for the purposes of adapting it to satisfy FFT needs?

Is there some other solution maybe?

We already discussed this a bit here.

[ckormanyos]

Is this crazy, but since we strive to use so many composite types in FFT, what if FFT were a part of Boost.Multiprecision? Turn it the other way around. Initially I do not like the feeling of that. I could also imagine writing, yet another, subset of something like complex solely for the FFT type(s) that are created within the inner parts of FFT.

[cosurgi]

It is not crazy. I was wondering myself if it should be a part of Multiprecision.

[cosurgi]

Math - by name - should be mathematical functions as I see it. Multiprecision - so far - are mostly numeric types with higher and higher precision.

Difficulties begin when you try to have higher and higher precision and mathematical functions. So far it worked, thanks to ADL, inside Boost.Math. Am I correct?

So if FFT is not in Boost.Multiprecision, but is in Boost.Math then heavy use of ADL and template specializations is a must.

What Boost users would say if FFT was in Multiprecision. I guess a surprise. We can rationalize this solution.

So which one is better. I seriously don't know.

[ckormanyos]

Is this crazy, but since we strive to use so many composite types in FFT, what if FFT were a part of Boost.Multiprecision?

Boost users would say if FFT was in Multiprecision. I guess a surprise. We can rationalize this solution.

If they had said and done nothing, they would have gotten it anyway. Sometime within the next 2 years or so I would have put some real-to-half-complex FFTs and their convolutions in a utilitarian part of Multiprecision on the road to millions of digits. In that particular case there would have been essentially no announcement from our side and all of a sudden it would be there in Boost.Multiprecision.

It is hard to think of another multiple-precision library that does not ship with its own home-brewed FFT.

@Lagrang3 does this make any problems go away? We7you could see how it might be possible to adapt the namespaces and arch for FFT to be within namespace boost::multiprecision if there is time?

[cosurgi]

I have managed to make_boost_complex for float128, cpp_bin_float and mpfr_float. It appears that cpp_dec_float and mpfi_float are missing some traits like signed_types, float_types, exponent_type. Or I did something wrong. Anyway here's my latest version of this code:

#include <boost/config.hpp>
#include <boost/cstdfloat.hpp>
#include <boost/type_traits/is_complex.hpp>
#include <boost/multiprecision/float128.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/multiprecision/mpfr.hpp>
#include <boost/multiprecision/mpfi.hpp>
#include <boost/core/demangle.hpp>
// now include complex types
#include <boost/multiprecision/complex128.hpp>
#include <boost/multiprecision/cpp_complex.hpp>
#include <boost/multiprecision/mpc.hpp>

namespace boost { namespace multiprecision {
template<typename T>
struct make_boost_complex {
    using type = std::complex<T>;
};

template<>
struct make_boost_complex<boost::multiprecision::float128> {
    using type = boost::multiprecision::complex128;
};
/*
template<>
struct make_boost_complex< number<backends::float128_backend, (boost::multiprecision::expression_template_option)0 > > {
    using type = boost::multiprecision::complex128;
};
*/
template <unsigned Digits, backends::digit_base_type DigitBase, class Allocator, class Exponent, Exponent MinExponent, Exponent MaxExponent, expression_template_option ExpressionTemplates>
struct make_boost_complex< number< cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>, ExpressionTemplates> > {
    using type = number<complex_adaptor<cpp_bin_float<Digits, DigitBase, Allocator, Exponent, MinExponent, MaxExponent>>, ExpressionTemplates>;
};

template <unsigned Digits, mpfr_allocation_type AllocationType, expression_template_option ExpressionTemplates>
struct make_boost_complex< number< backends::mpfr_float_backend<Digits, AllocationType>, ExpressionTemplates> > {
    using type = number<mpc_complex_backend<Digits>, ExpressionTemplates>;
};

template <typename T> struct is_boost_complex {
    static constexpr bool value = std::is_same<
                      typename make_boost_complex<typename T::value_type>::type
                    , typename std::decay<T>::type
                    >::value;
};
}} // boost::multiprecision

template<typename Real>
void print() {
    using std::sqrt;
    using Complex = typename boost::multiprecision::make_boost_complex<Real>::type;
    auto name_r = boost::core::demangle(typeid(Real).name());
    auto name_c = boost::core::demangle(typeid(Complex).name());
    std::cout << std::setw( 8) << sqrt(Complex(-1));
    std::cout << std::setw(12) << boost::is_complex<Complex>::value;
    std::cout << std::setw(18) << boost::multiprecision::is_boost_complex<Complex>::value << "  ";
    std::cout << name_c << std::endl;
//  std::cout << name_r << std::endl << std::endl;
}

template <typename... T>
constexpr inline auto for_each = [](auto&& f) {
    (f.template operator()<T>(), ...);
};

int main() {
    auto for_each_test_type = for_each<
              float
            , double
            , long double
            , boost::multiprecision::float128
            , boost::multiprecision::cpp_bin_float_50
            , boost::multiprecision::cpp_bin_float_quad
            , boost::multiprecision::cpp_dec_float_100
            , boost::multiprecision::number<boost::multiprecision::cpp_bin_float<200> >
            , boost::multiprecision::mpfr_float_100
            , boost::multiprecision::number<boost::multiprecision::mpfr_float_backend<200> >
            , boost::multiprecision::mpfi_float_100
            , boost::multiprecision::number<boost::multiprecision::mpfi_float_backend<200> >
        >;

    std::cout << "sqrt(-1) boost::is_complex boost::is_boost_complex type" << std::boolalpha << std::endl;
    for_each_test_type([]<typename T>() {
        print<T>();
    });

//  std::cout << boost::core::demangle(typeid(boost::multiprecision::cpp_dec_float_100).name()) << "\n";
//  std::cout << boost::core::demangle(typeid(boost::multiprecision::number<boost::multiprecision::complex_adaptor<boost::multiprecision::cpp_dec_float_100>,boost::multiprecision::et_off>).name()) << "\n";
//  std::cout << boost::core::demangle(typeid(boost::multiprecision::mpfi_float_100).name()) << "\n";
//  std::cout << boost::core::demangle(typeid(boost::multiprecision::number<boost::multiprecision::complex_adaptor<boost::multiprecision::mpfi_float_100>,boost::multiprecision::et_off>).name()) << "\n";
}

[cosurgi]

Forgot to mention that I compile with:

g++ test.cpp -o test -std=gnu++17 -lboost_system -lmpfr -lquadmath -Werror -Wformat -Wformat-security -Wformat-nonliteral -Wall -Wextra -Wnarrowing -Wreturn-type -Wuninitialized -Wfloat-conversion -Wcast-align -Wdisabled-optimization -Wtrampolines -Wpointer-arith -Wswitch-bool -Wwrite-strings -Wnon-virtual-dtor -Wreturn-local-addr -lmpfi -lmpc

[Lagrang3]

Now we need something which, in a universal way, constructs a Complex type for each of them. We have a hackish tool for that here already used in the code in several places.

@Lagrang3 already noted in the comments that it is needed.

We can actually get away without the complex selector. First of all the complex dft interface is a template on the complex type. The problem comes once we define the real dft interface with a real template parameter. If we want to convert the N real numbers to N complex numbers, there is no way we could know the type of the complex. I have two possible solution:

1. the rfft is parametrized on the complex type and not on the real type or
2. the rfft is parametrized on the real type, but any function that deals with complex numbers (such as real_to_complex and complex_to_real) is a template on the complex type.

[Lagrang3]

It is hard to think of another multiple-precision library that does not ship with its own home-brewed FFT.

I guess that, eventually, boost::multiprecision will have to choose between using boost::math::fft or duplicating it.

[ckormanyos]

eventually, boost::multiprecision will have to choose between using boost::math::fft or duplicating it

Indeed, which has, in fact, lead us to this issue.

[cosurgi]

eventually, boost::multiprecision will have to choose between using boost::math::fft or duplicating it

Indeed, which has, in fact, lead us to this issue.

And also issue #20, if Boost.Math and Boost.Multiprecision cannot depend on each other then to solve both this issue and #20 simultaneously maybe the solution is to create a separate library: Boost.FFT, then both Boost.Math and Boost.Multiprecision could depend on Boost.FFT in a clear manner.

1. the rfft is parametrized on the complex type and not on the real type or
2. the rfft is parametrized on the real type, but any function that deals with complex numbers (such as real_to_complex and complex_to_real) is a template on the complex type.

Indeed this would remove all dependencies altogether. A purely templated solution in Boost.FFT. I must admit that I like both of these. @ckormanyos what would you say?

Also about the 2nd approach: the complex type could be deduced from the type stored inside the target container.

[ckormanyos]

And also issue #20, if Boost.Math and Boost.Multiprecision cannot depend on each other then to solve both this issue and #20 simultaneously maybe the solution is to create a separate library: Boost.FFT, then both Boost.Math and Boost.Multiprecision could depend on Boost.FFT in a clear manner.

Yes that would have the advantage of clarity. I would be somewhat hesitant to tackle the overhead of creating/establishing/defending-the-need-for an entirely new library for the sole purpose of FFT transformations.

My strategy for Multiprecision is to make it also standalone. In that case, we could pull Multiprecision's FFT into Math and retain the standalone character.

If we put FFT in Multiprecision, a lot of clients might say, cool you finally have FFTs, but what does Multiprecision have to do with my engineering domain? That makes it kind of silly to have FFT in Multiprecision. That would bring us back to Math.

All this gets kind of circular, and you do ultimately end up at separate library thoughts. But I do still resist that.

My favorite at the state of our tech would be to get Multiprecision's dependencies low or standalone and include the tiny parts of Multiprecision in Math's interface/implementation for FFT, meaning FFT in math+make Multiprecision's dependencie(s) lean.

[cosurgi]

All this gets kind of circular, and you do ultimately end up at separate library thoughts. But I do still resist that.

I understand these reservations and fully support them.

My favorite at the state of our tech would be to get Multiprecision's dependencies low or standalone and include the tiny parts of Multiprecision in Math's interface/implementation for FFT, meaning FFT in math+make Multiprecision's dependencie(s) lean.

Yes, I think this is possible. Only small change needed in Boost.Multiprecision : integrate multiprecision_complex.hpp file. If I didn't miss anything, that should be totally enough. The only thing to decide is the name for boost::multiprecision::complex<…> template typename.

Unless by "tiny parts" you meant that this file should stay in Boost.Math, which is also fine for me.

Templatizing FFT solely on the Complex type is an alternative solution, which might make this file unnecessary. But somehow I still do like the notion of integrating all complex types available within Boost.Multiprecision into a single, easy to use, type.

[cosurgi]

I just noticed that maybe Boost.Math in general suffers this problem: "What complex type to use"? I just stumbled on an example with incorrect complex in return type in spherical harmonics. It forces std::complex ignoring complex_adaptor<cpp_bin_float<…>> or mpc_complex MPFR possibilities. It may be a good move to ask other Boost.Math authors what to do with this problem.

[Lagrang3]

I just noticed that maybe Boost.Math in general suffers this problem: "What complex type to use"? I just stumbled on an example with incorrect complex in return type in spherical harmonics. It forces std::complex ignoring complex_adaptor<cpp_bin_float<…>> or mpc_complex MPFR possibilities. It may be a good move to ask other Boost.Math authors what to do with this problem.

To me there is no such big step in moving from a real to complex. The complex is just two reals with more methods added. You see it was not difficult to implement a simple complex for the rfft. I wonder why the standard does not support std::complex for all arithmetic types. In the end, the compiler's implementation of std::complex (at least GCC) are very generic and they "eat" almost everything you throw in as T parameter.

[Lagrang3]

Containers do not ask many questions about the parameter T.

[ckormanyos]

maybe Boost.Math in general suffers this problem: "What complex type to use"?

wonder why the standard does not support std::complex for all arithmetic types.

Often times we have encountered this problem. At this point, I actually think we should consult some of the other Boost.Math developers. Where do you get complex for built-ins and UDTs?

Regarding the standard, I would perceive generic-complex as a potential improvement for the language. I kind of lost touch with the SGs and what's going on. I could ask around if anyone is working on this issue. I think SG6 (numerics) is the place to ask or elsewhere via proposal. I'm not sure at the moment if anyone is or not... But maybe if not, a proposal would be the right place to start moving forward on the standard.

[cosurgi]

@jzmaddock warned me that using std::complex on non fundamental floating (user defined types) point types is undefined behavior.

[cosurgi]

Link to other post where @ckormanyos mentions a bit controversial approach to this problem with cstdfloat_complex_std.hpp, which works with float_internal128_t but not complex_adaptor<cpp_bin_float<…>> and not with mpc_complex<…> MPFR type.

[cosurgi]

But maybe if not, a proposal would be the right place to start moving forward on the standard.

yes, definitely a standardized answer would do wonders here.

[ckormanyos]

mentions a bit controversial approach to this problem with cstdfloat_complex_std.hpp, which works with

Yes. I did mention that, and also mentioned its potential controversy and limitations. What I failed to mention is that the particular complex headers are not intended to be used directly. The intent is to simply #include <boost/cstdfloat.hpp>, ensure to compile with BOOST_MATH_USE_FLOAT128 and you can use Boost.Math's own version of boost::float128_t with Boost.Math's own version of std::complex<boost::float128_t>.

There are many drawbacks to this solution. And it does not handle higher multiprecision types such as those we need also like cpp-bin_float, etc.

The code below exhibits the abilities of <boost/cstdfloat.hpp>.

#include <complex>
#include <iomanip>
#include <iostream>

#include <boost/cstdfloat.hpp>

#if defined(BOOST_MATH_USE_FLOAT128)
using local_float_type = boost::float128_t;
#define MY_FLOAT_C(x) BOOST_FLOAT128_C(x)
#else
using local_float_type = boost::float64_t;
#define MY_FLOAT_C(x) BOOST_FLOAT64_C(x)
#endif

int main()
{
  local_float_type x(MY_FLOAT_C(1.1));
  local_float_type y(MY_FLOAT_C(2.2));

  std::complex<local_float_type> z(x, y);
  std::complex<local_float_type> s = sin(z);

  std::cout << std::setprecision(std::numeric_limits<local_float_type>::digits10) << s << std::endl;
}

Cc: @cosurgi and @Lagrang3

[ckormanyos]

But maybe if not, a proposal would be the right place to start moving forward on the standard.

yes, definitely a standardized answer would do wonders here.

I will query at SG6 and in the papers and try to find out if such a proposal is/was undertaken. This does not resolve the issue at hand "where do you get your complex", because standardization is typically a long-term process.

[cosurgi]

I just found a complex_result_from_scalar in:

1. number_base.hpp with comment: "// individual backends must specialize this trait." exactly our thoughts.
2. complex_adaptor.hpp
3. mpc.hpp

I discovered it accidentally in a call to boost::multiprecision::polar(…) which constructed a complex type. It wasn't used in spherical harmonics (mentioned above in which std::complex<UDT> was used.).

This complex_result_from_scalar is another contender to "where do you get your complex". I wonder if it could replace our boost::multiprecision::complex which we introduced in this branch.

[jzmaddock]

This complex_result_from_scalar is another contender to "where do you get your complex".

Within Multiprecision, yes please, it would be best to stick to all the same traits class for this.

[cosurgi]

Also we have a similar problem about recognizing whether a type boost::is_complex (it only recognizes std::complex) or is not complex. See https://github.com/BoostGSoC21/math/issues/12