Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Generic fft interface #8

Merged
merged 27 commits into from
Jul 5, 2021
Merged

Generic fft interface #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
104cce7
This needs a Fix: the #if (not defined(__clang__)) wrt. https://githu…
cosurgi Jun 23, 2021
2855074
different test tolerance for each backend
Jun 25, 2021
e0e3ee5
some comments
Jun 25, 2021
710ff43
generic DFT interface idea
Jun 25, 2021
7c22d19
Some L3 MSVC warns and use explicit logic ops
ckormanyos Jun 25, 2021
8c4a5f1
added complex root of unity functions
Jun 25, 2021
fd1f2b5
added tests for complex traits
Jun 25, 2021
6fc3a27
generic and complex DFT interface + tests
Jun 25, 2021
060b710
generic DFT: removed dependency on the neutral
Jun 25, 2021
72aa580
Merge branch 'my_generic_fft' into generic_fft_interface
Jun 25, 2021
24a56e2
fixed constexpr ctor compiler error
Jun 25, 2021
b9489a9
non-complex types supported + test of NTT
Jun 25, 2021
19ba639
Remove a few MSVC level 3 warnings and PR
ckormanyos Jun 26, 2021
d37da90
Merge branch 'develop' into generic_fft_interface
ckormanyos Jun 26, 2021
82faf77
Remove backend deps, warnings, spaces
ckormanyos Jun 26, 2021
fe927df
Try CI w/out BSL default but design still not right
ckormanyos Jun 26, 2021
ee44888
Re-try CI w/no default backend and MSVC errors
ckormanyos Jun 26, 2021
3fb3beb
first approach to convolutions
Jun 30, 2021
570d77d
convolution with padding
Jun 30, 2021
0bc8a27
Merge branch 'convolutions' into generic_fft_interface
Jun 30, 2021
5df6c0d
rader fft first implementation
Jul 1, 2021
392a2c4
fixed some bugs with rader and convolution impl.
Jul 1, 2021
5c2e26b
moved some discrete mathematics to dedicated header
Jul 3, 2021
8f4c355
power-mod using the template power function
Jul 3, 2021
4a33c19
naive prime factorization
Jul 3, 2021
69539f6
large primes to tests for accidental overflow
Jul 4, 2021
ac4b2dc
removed constexpr from couple of methods
Jul 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 101 additions & 28 deletions include/boost/math/fft.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,10 @@
#ifndef BOOST_MATH_FFT_HPP
#define BOOST_MATH_FFT_HPP

#include <algorithm>
#include <iterator>
//#include <boost/math/fft/fftw_backend.hpp>
//#include <boost/math/fft/gsl_backend.hpp>
#include <boost/math/fft/bsl_backend.hpp>
#include <vector>
#include <boost/math/fft/algorithms.hpp>

namespace boost { namespace math { namespace fft {

Expand All @@ -36,7 +36,7 @@
*/

template<typename RingType,
template<class U> class BackendType = bsl_dft >
template<class U> class BackendType>
// typename Allocator = Default_Allocator
class dft : public BackendType<RingType>
{
Expand Down Expand Up @@ -86,15 +86,15 @@
OutputIteratorType out,
typename std::enable_if<( (std::is_convertible<InputIteratorType, const RingType*>::value == true)
&& (std::is_convertible<OutputIteratorType, RingType*>::value == false))>::type* = nullptr)
{
{
resize(std::distance(in_first,in_last));
my_mem.resize(size());

if(ex==execution_type::backward)
backend_t::backward(in_first,my_mem.data());
else
backend_t::forward(in_first,my_mem.data());

std::copy(std::begin(my_mem), std::end(my_mem), out);
}

Expand All @@ -110,23 +110,26 @@
resize(std::distance(in_first,in_last));
my_mem.resize(size());
std::copy(in_first, in_last, std::begin(my_mem));

if(ex==execution_type::backward)
backend_t::backward(my_mem.data(),my_mem.data());
else
backend_t::forward(my_mem.data(),my_mem.data());

std::copy(std::begin(my_mem),std::end(my_mem), out);
}

public:
using backend_t = BackendType<RingType>;
using backend_t::size;
using backend_t::resize;

// complex types ctor. n: the size of the dft
constexpr dft(unsigned int n) : backend_t{ n } { }



// ring types ctor. n: the size of the dft, w: an n-root of unity
constexpr dft(unsigned int n, RingType w) : backend_t( n, w ) { }

template<typename InputIteratorType,
typename OutputIteratorType>
void forward(
Expand All @@ -135,7 +138,7 @@
{
execute(execution_type::forward,in_first,in_last,out);
}

template<typename InputIteratorType,
typename OutputIteratorType>
void backward(
Expand All @@ -147,43 +150,113 @@
};

// std::transform-like Fourier Transform API
template<template<class U> class backend = bsl_dft,
// for complex types
template<template<class U> class backend,
typename InputIterator,
typename OutputIterator>
void dft_forward(InputIterator input_begin,
InputIterator input_end,
OutputIterator output)
{
using input_value_type = typename std::iterator_traits<InputIterator >::value_type;
//using output_value_type = typename std::iterator_traits<OutputIterator>::value_type;

//static_assert(std::is_same<input_value_type, output_value_type>::value,
// "Input and output types mismatch");

dft<input_value_type, backend> plan(std::distance(input_begin, input_end));

dft<input_value_type, backend> plan(static_cast<unsigned int>(std::distance(input_begin, input_end)));
plan.forward(input_begin, input_end, output);
}

// std::transform-like Fourier Transform API
template<template<class U> class backend = bsl_dft,
// for complex types
template<template<class U> class backend,
class InputIterator,
class OutputIterator>
void dft_backward(InputIterator input_begin,
InputIterator input_end,
OutputIterator output)
{
using input_value_type = typename std::iterator_traits<InputIterator >::value_type;
//using output_value_type = typename std::iterator_traits<OutputIterator>::value_type;

//static_assert(std::is_same<input_value_type, output_value_type>::value,
// "Input and output types mismatch");

dft<input_value_type, backend> plan(std::distance(input_begin, input_end));

dft<input_value_type, backend> plan(static_cast<unsigned int>(std::distance(input_begin, input_end)));
plan.backward(input_begin, input_end, output);
}

// std::transform-like Fourier Transform API
// for Ring types
template<template<class U> class backend,
typename InputIterator,
typename OutputIterator,
typename value_type>
void dft_forward(InputIterator input_begin,
InputIterator input_end,
OutputIterator output,
value_type w)
{
using input_value_type = typename std::iterator_traits<InputIterator >::value_type;
dft<input_value_type, backend> plan(static_cast<unsigned int>(std::distance(input_begin, input_end)),w);
plan.forward(input_begin, input_end, output);
}

// std::transform-like Fourier Transform API
// for Ring types
template<template<class U> class backend,
class InputIterator,
class OutputIterator,
typename value_type>
void dft_backward(InputIterator input_begin,
InputIterator input_end,
OutputIterator output,
value_type w)
{
using input_value_type = typename std::iterator_traits<InputIterator >::value_type;
dft<input_value_type, backend> plan(static_cast<unsigned int>(std::distance(input_begin, input_end)),w);
plan.backward(input_begin, input_end, output);
}

template<template<class U> class backend,
typename InputIterator1,
typename InputIterator2,
typename OutputIterator>
void convolution(InputIterator1 input1_begin,
InputIterator1 input1_end,
InputIterator2 input2_begin,
OutputIterator output)
{
using input_value_type = typename std::iterator_traits<InputIterator1>::value_type;
using real_value_type = typename input_value_type::value_type;

const long N = std::distance(input1_begin,input1_end);
const long N_extended = detail::is_power2(N) ? N : detail::upper_bound_power2(2*N-1);

std::vector<input_value_type> In1(N_extended),In2(N_extended),Out(N_extended);

std::copy(input1_begin,input1_end,In1.begin());

InputIterator2 input2_end{input2_begin};
std::advance(input2_end,N);
std::copy(input2_begin,input2_end,In2.begin());

// padding
for(long i=N;i<N_extended;++i)
In1[i]=In2[i]=input_value_type{0};

// fake N-periodicity
if(N!=N_extended)
for(long i=1;i<N;++i)
In2[N_extended-N+i] = In2[i];

dft<input_value_type, backend> plan(static_cast<unsigned int>(N_extended));
plan.forward(In1.begin(),In1.end(),In1.begin());
plan.forward(In2.begin(),In2.end(),In2.begin());

// direct convolution
std::transform(In1.begin(),In1.end(),In2.begin(),Out.begin(),std::multiplies<input_value_type>());

plan.backward(Out.begin(),Out.end(),Out.begin());

const real_value_type inv_N = real_value_type{1}/N_extended;
for(auto & x : Out)
x *= inv_N;

std::copy(Out.begin(),Out.begin() + N,output);
}

} } } // namespace boost::math::fft

#endif // BOOST_MATH_FFT_HPP
22 changes: 22 additions & 0 deletions include/boost/math/fft/abstract_ring.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@
#include <boost/multiprecision/complex128.hpp>
#endif
#include <boost/multiprecision/complex_adaptor.hpp>
#include <boost/math/constants/constants.hpp>

/*
RingType axioms:
Expand All @@ -50,6 +51,27 @@
*/

namespace boost { namespace math { namespace fft { namespace detail {

#if __cplusplus < 201700L
template<typename ...> using void_t = void;
#else
using std::void_t;
#endif


template<class T, typename = void_t<> >
struct is_complex : std::false_type
{};

template<class T >
struct is_complex<T,
void_t<
typename T::value_type,
decltype(sin(std::declval<typename T::value_type>())),
decltype(cos(std::declval<typename T::value_type>()))
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm.. if I get the intentions correctly it should be checking for presence of

decltype(sin(std::declval<T>())),
decltype(cos(std::declval<T>()))

Or did I misunderstood something?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I see, you are computing sin(…) and cos(…) of phase ∈ ℝ. OK. I will need to ponder this for a while.

Copy link
Collaborator

@Lagrang3 Lagrang3 Jun 28, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, it is not so important to compute the sin/cos of the complex number. Our implementation needs the sin/cos of the underlying real phase.

> > : std::true_type
Copy link
Collaborator

@cosurgi cosurgi Jun 27, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's a very nice template. However I am a little bit concerned with it's name. Maybe better to be more explicit that it checks for presence of sin(…) and cos(…). Like maybe has_sin_cos<…> ?

Following this, I think we should to go through Boost libs and add tests for stuff which happens to have sin(…) and cos(…), this one is the first contender.

It's also possible to calculate the sin(…) and cos(…) of a matrix.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So far we have two kinds of template algorithms:

  1. those that assume the type is complex so that we can use sin/cos to compute the roots of unity (so that we don't have to ask for the user for the roots and we don't have to use power functions to derive smaller roots from larger ones);
  2. and those who would work even in a larger domain (this domain include complex as well) of algebraic ring types.

The intention of this template is to tell if the argument represents a complex number, so that we can choose one DFT implementation or the other. Checking for sin/cos is a first approach to deduce the answer for unknown types. We could change the logic if needed to.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, I see. Then maybe simpler it would be to use the native boost::is_complex from #include <boost/type_traits/is_complex.hpp> ?

Copy link
Collaborator

@Lagrang3 Lagrang3 Jun 28, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That would be great! But the documentation says:
"If T is a complex number type then true (of type std::complex for some type U), otherwise false."
That means (I've just tested)
boost::is_complex<boost::multiprecision::cpp_complex_quad>::value
is false.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure, we could alias and specialize the template boost::is_complex inside the boost::math::fft namespace to evaluate to true all boost::multiprecision::cpp_complex<...> types. But then I guess boost::math would be dependent on boost::type_traits and boost::multiprecision.

Copy link
Collaborator

@cosurgi cosurgi Jun 28, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That means (I've just tested) boost::is_complex<boost::multiprecision::cpp_complex_quad>::value is false.

@ckormanyos , this seems like a bug to me, what would you say about this?

EDIT: or maybe there is some other is_complex template for checking multiprecision stuff, like multiprecision::is_complex?

Then boost::is_complex<…> or boost::multiprecision::is_complex<…> would do the trick :)

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As for the Type Traits library, that's not a bug
https://github.com/boostorg/type_traits/blob/develop/include/boost/type_traits/is_complex.hpp
The value is true if the type is std::complex<...> otherwise false. Hence the template tells if the type is "complex", where "complex" means std::complex.

Multiprecision has its own definition of the is_complex template
https://github.com/boostorg/multiprecision/blob/develop/include/boost/multiprecision/traits/is_complex.hpp
But it does the same thing.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for checking this. We will have time figure this out. Maybe for example check if the type uses complex_adaptor<…>. But it's not urgent. What matters is that now I understand your intentions in this part of code :)

Copy link
Member Author

@ckormanyos ckormanyos Jun 29, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In complex_adaptor, we use things like

enable or test a Boolean value from:

boost::multiprecision::detail::is_complex<Result>::value

I need to check again if this works for non-built-in (i.e., multiprecision) types.

But here is a rather restricting constraint. Very much work was done to make Boost.Math standalone, requiring no other Boost library dependencies other than itself. It might be incongruent with the evolution of Boost.Math to have FFT depend on Boost.Multiprecision.

Copy link
Collaborator

@cosurgi cosurgi Jun 29, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

.. to have FFT depend on Boost.Multiprecision.

I definitely agree with that. We will need to do more SFINAE (or #ifdef :) ) tricks to make sure that a simple line that calls FFT on a simple double type does not pull extra dependencies with it.

{};


// Recognizes:
// → fundamental types
Expand Down
Loading