api.tex

\section{Detailed API}
Boost.Math.FFT brings two kinds of interfaces: one-time and planned transforms.
The first one is inspired by the standard library's \texttt{std::transform} 
and the second by \fftw's plan interface. Listings
\ref{ls:complex_dft_transform} and \ref{ls:complex_dft_plan} are two examples of
the use of the one-time and planned interfaces respectively.

The planned transform is used with
the declaration of a plan variable and subsequently calling its methods.
The internal state of the plan is determined by the size of the \dft\ transform,
and throughout its life it caches the necessary data and memory allocation
needed to perform \dft s on request. In principle it should be faster to call a
method of a plan, than the equivalent routine from the one-time interface.

The library's interface can be further classified according to the domain where
the \dft\ is applied: let it be complex numbers, real and halfcomplex, or a
purely algebraic set satisfying the ring axioms. This distinction is necessary
because all three behave differently and the user expects different interfaces.

Consider the complex case. There are only two types of \dft\ we can perform on
a $\Complex^n$ vector, they are the \verb|forward| and the \verb|backward|, 
because there are only two roots of unity in that domain:
$w = \exp( - 2 \pi i/n)$ 
and
$w^{-1} = \exp( 2 \pi i/n)$.
While in general for an algebraic transform on a type that satisfies the ring
axioms the \dft\ engine has no way of knowing \emph{a priori} what is the root
of unity $w$ and that one must be supplied by the user either to the constructor
of the plan or as an extra argument to the one-time transform routine. In this
case, once $w$ is known there are again only two types of \dft s: \verb|forward| 
when we use $w$ or \verb|backward| when using $w^{-1}$.

The real/halfcomplex interface works very differently from the complex and the
algebraic. Since the result of any \dft\ on a $\Real^n$ vector yields always a
vector in $\hComplex^n$ and viceversa; this interface contains the methods:
\verb|real_to_halfcomplex| and \verb|halfcomplex_to_real| to emphasize that
property. Because of the trivial relations between the results of the forward
and backward transforms applied to real and/or halfcomplex data (see fig.
\ref{fig:halfcomplex}) we do not provide a choice for forward or backward
transform. \verb|real_to_halfcomplex| will in effect compute a forward \dft\ on
a $\Real^n$ vector and \verb|halfcomplex_to_real| will compute a backward \dft\ on a
$\hComplex^n$ vector. We also provide the methods \verb|real_to_complex|
and \verb|complex_to_real| similar to the previous \verb|real_to_halfcomplex|
and \verb|halfcomplex_to_real| respectively, where the halfcomplex part is
represented in an array of $n$ complex numbers, hence keeping redundant values.

The library computes the requested operations using one of three backends:
\gsl, \fftw\ or our own Boost licensed engine that we call \bsl. 

\subsection{Complex \dft}
Complex valued \dft\ are best handled by Boost.Math.FFT complex oriented
interface. All three backends are accessible through the plan-like classes:
\verb|bsl_dft|, \verb|gsl_dft| and \verb|fftw_dft|, representing the three
engines \bsl, \gsl\ and \fftw\ respectively. They share the same interface hence
it is enough to describe only the member functions of \verb|bsl_dft|.
\begin{lstlisting}[language=C++,caption=Complex plan declaration.]
template< class Complex = std::complex<double>, 
          class Allocator =  std::allocator<Complex> >
class bsl_dft;
\end{lstlisting}

Notice that even though the complex \dft\ plans are template on the complex
type, not all backends handle all complex types. \gsl\ only works for
\verb|std::complex<double>| and \fftw\ work with a richer collection of
precision types
\verb|std::complex<float>|, \verb|std::complex<double>|, 
\verb|std::complex<long double>| and
\verb|boost::multiprecision::complex128|\footnote{or equivalently
\texttt{boost::multiprecision::complex<boost::multiprecision::float128>}.}.
On the other hand, the implementation of \verb|bsl_dft| allows, in principle, to
compute \dft s
for any complex type as long as it has the essential complex
algebraic methods plus sine and cosine functions for its internal real
\verb|value_type|.
\verb|bsl_dft| can handle, for instance besides the previously mentioned complex types,
also the family of \verb|boost::multiprecision|'s \verb|cpp_complex|
and \verb|mpc_complex| based on \verb|cpp_bin_float| and
\verb|mpfr_float_backend| real types respectively; eg.
\verb|cpp_complex_50|, 
\verb|cpp_complex_quad| and 
\verb|mpc_complex_50|.

\subsubsection*{Member functions}
\begin{lstlisting}[language=C++,caption=Complex plan constructor.,label=ls:complex_plan]
template< class Complex, class Allocator >
constexpr bsl_dft<Complex,Allocator>::bsl_dft(
            unsigned int n, const Allocator& alloc = Allocator{} );
\end{lstlisting}
The constructor of a complex \dft\ plan takes two arguments (listing
\ref{ls:complex_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).
\begin{lstlisting}[language=C++,caption=Complex plan forward function.]
template<class InputIt,
         class OutputIt>
void forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex plan backward function.]
template<class InputIt,
         class OutputIt>
void backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex plan size.]
constexpr std::size_t size() const;
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Resize.]
void resize(std::size_t new_size);
\end{lstlisting}

\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.
\begin{lstlisting}[language=C++,caption=Forward complex transform.]
template<class InputIt,
         class OutputIt>
void bsl_transform::forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Backward complex transform.]
template<class InputIt,
         class OutputIt>
void bsl_transform::backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}

\subsubsection*{Examples}
See listings \ref{ls:complex_dft_transform} and \ref{ls:complex_dft_plan} 
for examples of the use of the complex one-time transform and plan interfaces.

\subsection{Algebraic \dft}
Algebraic type \dft\ is understood in the sense of definition \ref{def:dft}.
It can be applied to algebraic types that satisfy the ring axioms (see appendix
\ref{app:ring}) as long as a root of unity is known for the required \dft\ size.
Boost.Math.FFT has achieved this level of abtraction with the use of C++
templates. There is not equivalent functionality in \gsl\ or \fftw\ libraries
for this kind of transform, hence we only describe the \bsl\ plan-like backend
called \verb|bsl_algebraic_dft|.
\begin{lstlisting}[language=C++,caption=Algebraic plan.,label=ls:alg_plan]
template< class T = std::complex<double>, 
          class Allocator = std::allocator<T> >
class bsl_algebraic_dft;
\end{lstlisting}
The constructor of an algebraic \dft\ plan takes three arguments (listing
\ref{ls:alg_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is an $N$-root of unity $w$ and the third is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).

\subsubsection*{Member functions}
\begin{lstlisting}[language=C++,caption=Algebraic plan constructor.]
template< class T, class Allocator >
constexpr bsl_algebraic_dft<T,Allocator>::bsl_algebraic_dft(
            unsigned int n, T w, const Allocator& alloc = Allocator{} );
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Algebraic plan forward function.]
template<class InputIt,
         class OutputIt>
void forward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Algebraic plan backward function.]
template<class InputIt,
         class OutputIt>
void backward(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Size.]
constexpr std::size_t size() const;
\end{lstlisting}

% The fact is that the algebraic_dft does have a resize method, but resizing the
% DFT here must be followed with a change of the internal root of unity w. A
% correct resize method here requires two arguments: the new size and the new
% value of w.
% \begin{lstlisting}[language=C++]
% void resize(std::size_t new_size);
% \end{lstlisting}

\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.

\begin{lstlisting}[language=C++,caption=Forward algebraic transform.]
template<class InputIt,
         class OutputIt,
         class T>
void bsl_algebraic_transform::forward(
       InputIt first1, InputIt last1, OutputIt d_first, const T w);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Backward algebraic transform.]
template<class InputIt,
         class OutputIt,
         class T>
void bsl_algebraic_transform::backward(
       InputIt first1, InputIt last1, OutputIt d_first, const T w);
\end{lstlisting}

\subsubsection*{Examples}
Complex numbers are a particular case of an algebraic ring, which fortunately we
know for every $N$ only two possible $N$-roots of unity; ie. 
$\exp(-i \frac{2\pi}{N})$
and
$\exp(i \frac{2\pi}{N})$, where $i=\sqrt{-1}$. Hence we could use the algebraic
interface to compute complex \dft s. Listing \ref{ls:complex_alg_dft} shows how
could we compute the \dft\ of a complex array using the algebraic interface.
\begin{lstlisting}[language=C++,caption=Complex DFT using the generic interface.,label=ls:complex_alg_dft]
using Real = double;
using Complex = std::complex<Real>;

std::vector<Complex> A{1.,-1.,3.,.5,9.},B,C;
const int N = A.size();
const Real w_phase = 2*boost::math::constants::pi<Real>()/N;
const Complex w{cos(w_phase), -sin(w_phase)}; // w is an N-root of unity in the complex

// compute the complex forward DFT using the algebraic interface
bsl_algebraic_transform::forward(A.cbegin(),A.cend(),std::back_inserter(B), w);
// compute the complex forward DFT using the complex interface
bsl_transform::forward(A.cbegin(),A.cend(),std::back_inserter(C));
// B and C are equal up to floating point errors
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Integer fast multiplication.]
using M_int = fft::my_modulo_lib::mint<Z337>;
// 85 is a primitive root of 337,
// ie. the smallest number k such that 85^k mod 337 = 1 is k=336
const M_int w{85};
const M_int inv_8{M_int{8}.inverse()};

// Multiplying 1234 times 5678 = 7006652
std::vector<M_int> A{4, 3, 2, 1, 0, 0, 0, 0};
std::vector<M_int> B{8, 7, 6, 5, 0, 0, 0, 0};

// forward FFT
boost::math::fft::bsl_algebraic_dft<M_int> P(A.size(),w);

P.forward(A.cbegin(),A.cend(),A.begin());
P.forward(B.cbegin(),B.cend(),B.begin());

// convolution in Fourier space
std::vector<M_int> AB;
std::transform(A.begin(), A.end(), B.begin(),
               std::back_inserter(AB),
               [](M_int x, M_int y) { return x * y; });

// backwards FFT
P.backward(AB.cbegin(),AB.cend(),AB.begin());
std::transform(AB.begin(), AB.end(), AB.begin(),
               [&inv_8](M_int x) { return x * inv_8; });

// carry the remainders in base 10
std::vector<int> C;
M_int r{0};
for (auto x : AB)
{
  auto y = x + r;
  C.emplace_back(int(y) % 10);
  r = M_int(int(y) / 10);
}
// C[] = 7006652
\end{lstlisting}

\subsection{Real \dft}
Real valued \dft\ are best handled by Boost.Math.FFT real oriented
interface. The plan-like classes specialized for real numbers are denoted
\verb|bsl_rdft|, \verb|gsl_rdft| and \verb|fftw_rdft|, 
representing the three
engines \bsl, \gsl\ and \fftw\ respectively. They share the same interface hence
it is enough to describe only the member functions of \verb|bsl_rdft|.
The real methods take advantage of the symmetries involved in the real valued
\dft s described in Subsection \ref{subsec:real}. In principle one could compute
those real transforms using complex numbers and the complex interface (eg. by
setting to zero the imaginary part of the input and ignoring the halfcomplex
symmetry), but the real routines will perform
half of the operations and occupy half the space in memory than the complex
unconstrained counterpart. 
\begin{lstlisting}[language=C++,caption=Real plan.]
template< class Real = double, 
          class Allocator =  std::allocator<Real> >
class bsl_rdft;
\end{lstlisting}
Real plans, as well as the complex plans seen previously, are templates on the
type that represents the real numbers. Again \gsl\ backend will work only for
\verb|double| representations while \fftw\ supports \verb|float|, \verb|double|,
and \verb|long double|, and \verb|boost::multiprecision::float128|. 
The \bsl\ functions are all templates on the type,
hence it will accept any type representing real numbers.

\subsubsection*{Member functions}

\begin{lstlisting}[language=C++,caption=Real plan constructor.,label=ls:real_plan]
template< class Real, class Allocator >
constexpr bsl_rdft<Real,Allocator>::bsl_rdft(
            unsigned int n, const Allocator& alloc = Allocator{} );
\end{lstlisting}
The constructor of a real \dft\ plan takes two arguments (listing
\ref{ls:real_plan}): the first is an integer number $N$ that specifies the size of
the transform, the second is a
reference to an allocator (allocators are discussed in section
\ref{subsec:allocators}).
\begin{lstlisting}[language=C++,caption=Real to halfcomplex.]
template<class InputIt,
         class OutputIt>
void real_to_halfcomplex(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Real to complex.]
template<class InputIt,
         class OutputIt>
void real_to_complex(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Halfcomplex to real.]
template<class InputIt,
         class OutputIt>
void halfcomplex_to_real(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Complex to real.]
template<class InputIt,
         class OutputIt>
void complex_to_real(InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Size.]
constexpr std::size_t size() const;
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=Resize.]
void resize(std::size_t new_size);
\end{lstlisting}

\subsubsection*{One-time transform methods}
The same functionality is accessible through the one-time transform methods
below.

\begin{lstlisting}[language=C++,caption=One-time transform real to halfcomplex.]
template<class InputIt,
         class OutputIt>
void bsl_real_transform::real_to_halfcomplex(
       InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform real to complex.]
template<class InputIt,
         class OutputIt>
void bsl_real_transform::real_to_complex(
       InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform halfcomplex to real.]
template<class InputIt,
         class OutputIt>
void bsl_real_transform::halfcomplex_to_real(
       InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}
\begin{lstlisting}[language=C++,caption=One-time transform complex to real.]
template<class InputIt,
         class OutputIt>
void bsl_real_transform::complex_to_real(
       InputIt first1, InputIt last1, OutputIt d_first);
\end{lstlisting}

\subsubsection*{Examples}
\begin{lstlisting}[language=C++,caption=Polynomial multiplication.]
using Real = double;
std::vector<Real> A{1.,4.,-5.,1.,0.,0.,0.,0.};
std::vector<Real> B{-1.,1.,2.,3.,0.,0.,0.,0.};
const std::size_t N = A.size();
std::vector< std::complex<Real> > TA(N),TB(N);

boost::math::fft::bsl_rdft<Real> P(N); 
P.real_to_complex(A.begin(),A.end(),TA.begin());
P.real_to_complex(B.begin(),B.end(),TB.begin());

for(unsigned int i=0;i<N;++i)
  TA[i]*=TB[i];
  
std::vector<Real> C(N);
P.complex_to_real(TA.begin(),TA.end(),C.begin());
std::transform(C.begin(), C.end(), C.begin(),
               [N](Real x) { return x / N; });
// C[] = [-1, -3, 11, 5, 3, -13, 3, 0,]
\end{lstlisting}

\subsection{Convolution}
Currently, there is no convolution transform implemented as a member of any of
the plan-like classes. Instead the circular convolution (see definition
\ref{def:convolution}) of two series of complex
numbers is handled by the one-time transform described in listing~\ref{ls:convolution}.
    
\begin{lstlisting}[language=C++,caption=Convolution.,label=ls:convolution]
template<typename InputIt1,
         typename InputIt2,
         typename OutputIt>
void bsl_transform::convolution(
       InputIt1 first1, InputIt1 last1, InputIt2 first2, OutputIt d_first);
\end{lstlisting}

\subsection{The one-time transform generalities}
For each plan type and member function there is an equivalent one-time transform
method. These are static methods of a template class denoted \verb|transform|. 
This design choice allows to select one backends programmatically.
For instance \verb|bsl_transform| is an alias for \verb|transform<bsl_dft<>>|,
\verb|bsl_real_transform| is and alias for \verb|transform<bsl_rdft<>>|, etc.
Table \ref{tab:transform} summarizes the naming of the one-time transform family.
\begin{table}[h]
    \centering
    \begin{tabular}{ll}
        Alias & Instantiation \\
        \hline
        \verb|bsl_transform| & \verb|transform< bsl_dft<> >|\\
        \verb|bsl_algebraic_transform| & \verb|transform< bsl_algebraic_dft<> >|\\
        \verb|bsl_real_transform| & \verb|transform< bsl_rdft<> >|\\
        \verb|gsl_transform| & \verb|transform< gsl_dft<> >|\\
        \verb|gsl_real_transform| & \verb|transform< gsl_rdft<> >|\\
        \verb|fftw_transform| & \verb|transform< fftw_dft<> >|\\
        \verb|fftw_real_transform| & \verb|transform< fftw_rdft<> >|\\
    \end{tabular}
    \caption{Summary of the one-time transform classes.}
    \label{tab:transform}
\end{table}

\subsection{The allocator template parameter}\label{subsec:allocators}
The plan-like interface has a template parameter named \verb|Allocator|.
The plan uses the \emph{allocator} supplied through the constructor, to allocate
all the internal heap memory necessary to compute \dft s.
Notice that very often, \dft s cannot be performed in place, and scratch space
is used. Furthermore, the work space is required to be a contiguous array in
memory and there is no way to deduce that property for any iterator before
C++20's \emph{contiguous iterator} concept; hence scratch space is used even in
cases when the in-place \dft\ could be possible in principle.

The example in listing \ref{ls:pmr_allocator}, shows how could the allocator
parameter could be used to harness memory allocation.

\begin{lstlisting}[language=C++,caption=Example use of custom allocators.,label=ls:pmr_allocator]
const std::size_t N = 100;
std::array<char,30000> buf;
boost::container::pmr::monotonic_buffer_resource
  pool{buf.data(),buf.size(),boost::container::pmr::null_memory_resource()};

using Real = double;
using Complex = boost::multiprecision::complex<Real>;
using Allocator = boost::container::pmr::polymorphic_allocator<Complex>;

std::vector<Complex,Allocator> A(N,Complex(),&pool);
std::vector<Complex,Allocator> B(N,Complex(),&pool);
// ...
bsl_dft<Complex,Allocator> plan(N,&pool);
plan.forward(A.begin(),A.end(),B.begin());
\end{lstlisting}